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Abstract 

We present a review of the current state of the art of cosmological dark matter simulations, 
with particular emphasis on the implications for dark matter detection efforts and studies of dark 
energy. This review is intended both for particle physicists, who may find the cosmological 
simulation literature opaque or confusing, and for astro-physicists, who may not be familiar with 
the role of simulations for observational and experimental probes of dark matter and dark energy. 
Our work is complementary to the contribution by M. Baldi in this issue, which focuses on the 
treatment of dark energy and cosmic acceleration in dedicated N-body simulations. 

Truly massive dark matter-only simulations are being conducted on national supercomputing 
centers, employing from several billion to over half a trillion particles to simulate the formation 
and evolution of cosmologically representative volumes (cosmic scale) or to zoom in on indi- 
vidual halos (cluster and galactic scale). These simulations cost millions of core-hours, require 
tens to hundreds of terabytes of memory, and use up to petabytes of disk storage. Predictions 
from such simulations touch on almost every aspect of dark matter and dark energy studies, and 
we give a comprehensive overview of this connection. We also discuss the limitations of the 
cold and collisionless DM-only approach, and describe in some detail efforts to include differ- 
ent particle physics as well as baryonic physics in cosmological galaxy formation simulations, 
including a discussion of recent results highlighting how the distribution of dark matter in halos 
may be altered. We end with an outlook for the next decade, presenting our view of how the field 
can be expected to progress. 
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1. Introduction 

For almost 80 years now [ 1| astronomers have been accumulating evidence that the dominant 
form of matter in the universe is dark and non-baryonic. Just a bit more than a decade ago, 
we obtained solid observational evidence of cosmic acceleration [3), requiring yet another 
mysterious contribution to the total energy budget of the universe, dark energy. The conclusion 



Email addresses: mqkSastro . berkeley . edu (Michael Kuhlen), mvogelsbScf a . harvard . edu (Mark 
Vogelsberger), reanguloOslac . Stanford, edu (Raul Angulo) 

Preprint submitted to Physics of the Dark Universe October 29, 2012 



is staggering: we live in a universe that is energetically dominated by dark matter (DM) and 
dark energy (DE). For DM at least we have a well-motivated theoretical framework, in which 
it is comprised of fundamental particles predicted to arise in extensions of the standard model 
of particle physics like supersymmetry [4|. There are many ongoing and proposed efforts to 
obtain experimental confirmation of the hypothesis of particle DM, for example through so called 
indirect (annihilation/decay) and direct (nuclear scattering) detection. These signatures depend 
on the detailed distribution of DM throughout the universe, from cosmic all the way down to 
sub-galactic scales. For DE, on the other hand, we have only a very rudimentary theoretical 
understanding, and are a firmly in an exploratory phase, conducting and designing future surveys 
and measurements that will provide us with additional clues to its nature. 

The last 40 years have seen tremendous progress in our understanding of cosmic structure and 
galaxy formation. Much of this understanding has come at the hand of beautifully simple analytic 
arguments and insights. The calculation of the Cold Dark Matter (CDM) power spectrum J5] [6], 
Press & Schechter theory [7 |, the statistics of peaks in Gaussian random fields [8], and White 
& Rees' galaxy formation model [9| are a few seminal examples. Yet, it is clear that purely 
analytical approaches have arrived at the limits of their reach. Fueled by continuing advances in 
numerical methods and computational capabilities, the future of structure formation and galaxy 
formation theory is going to be led by numerical simulations. 

Indeed, in the last few decades Moore's law combined with heavy infrastructure investments 
has already tremendously increased available computing resources, and the field of computa- 
tional cosmology has taken full advantage. Full-box simulations of substantial fractions of the 
observable Universe are now being conducted with over half a trillion particles, and zoom-in 
simulations of individual halos have exceeded the billion particle level. In this review we present 
a snapshot of the current state of the art of cosmological numerical simulations, with a particular 
emphasis on the DM and DE problems. Its intended audience is both the particle physicist with 
an interest in DM and DE, who may find the simulation literature to be opaque and confusing, 
as well as the astro-physicist, who may not be up to speed with observational and experimental 
probes of DM and DE or with the importance of simulations for their interpretation. Our work is 
complementary to the contribution by M. Baldi in this issue [10], which focuses on the treatment 
of DE and cosmic acceleration in dedicated N-body simulations. 

This review is structured as follows. We begin in ^2] with a review of the domain of DM 
simulations. What sorts of predictions do they make? How are these predictions relevant for 
DM detection efforts and for probes of DE? We cover astrophysical probes, indirect, and direct 
detection of DM, as well as probes of DE in the form of baryon acoustic oscillations, redshift- 
space distortions, cluster mass functions, and weak gravitational lensing. In ^3] we present a 
survey of the state of the art in late 2012 of cosmological DM simulations on cosmic, cluster, 
and galactic scales. We give an overview of the most commonly employed numerical techniques, 
then describe some of the largest simulations run to date, and the computational resources they 
used. We discuss the limits of the cold and collisionless DM-only approach, and efforts to go 
beyond it by simulating alternative DM physics. We include an discussion of hydrodynamic 
galaxy formation simulations, cover numerical methods, algorithmic and technical difficulties, 
and highlight some recent results regarding how the DM distribution in halos may be altered by 
baryonic physics. Finally, in |4]we present our vision of this field for the next decade. What are 
the important questions to tackle, and how best to do so? What developments should be pursued 
in order to take advantage of technological advances? 
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Figure 1: A 2 (k) a An(kl2n) i P(k), the linear power spectrum of density fluctuations at z = 0. The solid line is the 
canonical cold DM model with an Eisenstein & Hu (1997) 1 1 Q transfer function. The dashed line is a thermal relic warm 
DM model with mwDM = 8 keV |12| . The dotted line is an atomic DM model [13]. We used WMAP7 cosmological 
parameters O, f2„, = 0.265, H A = 0.735, Cl b = 0.0449, h = 0.71, cr g = 0.801, and n s = 0.963. 



2. Dark Matter Simulations and the Dark Universe 

The numerical simulation discussed in this review together span an enormous range of length 
scales, more than 8 orders of magnitude reaching from near horizon scale (~ 20 Gpc) down to 
sub-Galactic (tens of pc). Individually they focus on different regimes (see §|3]and Table[2]), but 
all have in common that they evolve the growth of DM density fluctuations all the way to the 
present epoch at redshift zeroQ 

The shape of the CDM power spectrum results in a hierarchical, bottom-up process of struc- 
ture formation, in which small and low mass objects collapse first and over time merge to form 
ever more massive structures, until the onset at z » 1 of DE induced accelerated expansion begins 
to halt further collapse. In Fig. [T] we show a plot of the linear dimensionless matter power spec- 
trum A 2 (k) = 47r(k/2nfP(k) at z = versus the wavenumber k of the fluctuation. Where A > 1, 
gravitational collapse will have proceeded to the non-linear regime and typical objects of the cor- 
responding mass will have collapsed. Cosmic scales, including the Baryon Acoustic Oscillation 



feature discussed in § 2.3 i, remain in the linear or mildly non-linear regime, while cluster and 
galactic scales are strongly non-linear. Note that computational demands grow strongly with the 
degree of non-linearity resolved in the simulation. 

Observational constraints from the Cosmic Microwave Background, cluster abundances, galaxy 
clustering, weak lensing and the Lyman-or forest have constrained the power spectrum of density 



'We deliberately omit from our discussion multi-billion particle simulations that focus only on the first billion years 
of cosmic evolution, for studying the epoch of reionization 1151 or early supermassive black hole growth 1161 . 
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fluctuations and provide a remarkably good agreement with the predictions of ACDM cosmol- 
ogy. On smaller scales (shaded gray in Fig. [TJ we currently do not have robust observational 
constraints, and here numerical simulations typically rely on extrapolations under the assump- 



tion of CDM theory. As we discuss in 4 3.3 different assumptions are plausible and are the 
subject of many ongoing investigations. As an example, we show two alternative models with a 
suppression of small scale power, a warm DM lfl2l and an atomic DM lfl3l model. 

2.1. Domain 

In the following we provide a brief summary of the domain of cosmological DM simulations, 
roughly organized from large scales to smalljj This is not meant to be an exhaustive review of 
all current results, but rather an overview of those results with particular relevance for DM and 
DE experiments. The references that we list provide a jumping-off point for further reading. 

i) Large Scale Structure 

The largest scale density fluctuations in the universe never grow beyond mildly non-linear, 
and even early CDM simulations predicted that the large scale distribution of DM in the 
universe is not completely homogeneous, instead exhibiting voids, walls, and filaments 
whose statistical description is in remarkable agreement with the large scale distribution 
of galaxies [e.g. [18]. 

Simulations spanning an ever larger fraction of the volume of the observable Universe at 
increasingly high resolution have been able to quantify the DM density and velocity fields 
as well as the halo mass function together with the full hierarchy of halo correlation 
functions, and their evolution with cosmic time |[T9l - l23"l . 

ii) Individual isolated halo properties 

On the scale of individual halos, DM-only numerical simulations have measured halo 
shapes to show significant departures from sphericality, with halos typically being prolate 
and increasingly so towards their centers. Major- to-minor axis ratios of 2 or greater are 
not uncommon, and more massive halos tend to be less spherical than lower mass halos 
11241 l25ll . Shapes and kinematics seem to be closely connected. While the spherically 
averaged anisotropy profile (J3(r) = 1 - 0.5 crf/o~$) grows from zero (isotropic) to about 
0.4 (mild radial anisotropy) |26 27 1, the local f3 values correlate with halo shape: positive 
(radial) on the major axis and negative (tangential) on the minor axis Il28l[29l . 

The DM mass distribution within halos is well described by a near-universal density pro- 
file, the so-called NFW profile |30|, which has the form of a double-power-law with the 
logarithmic slope y = dlogp/dlogr transitioning at the scale radius r s from y = -3 at 
large radii to y — - 1 in the center. More recent higher resolution simulations, however, 
have found a central slope shallower than y = — 1, indicating that the density profile may 
be better described by a functional form with a central slope gradually flattening to y — 0, 
e.g. the Einasto profile OTl l32ll . 

The scaling of the transition radius r s with halo mass, formation time, and environment is 
typically described in terms of a "concentration", defined as the ratio of the virial radius 
to the scale radius, c - Rvn/r s . DM simulations have quantified the concentration-mass 



2 Also see 1171 for a recent review of simulation results concerning the evolution and structure of CDM halos. 
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relationship, its scatter, and its evolution with time l33H35ll . Concentrations typically 
increase for lower mass halos, presumably reflecting their earlier collapse times when the 
mean density of the universe was higher, although recent work has reported an upturn of 
concentrations at high masses [36 1 presumably caused by out-of-equilibrium systems l37l . 

Lastly, we mention the remarkable finding from simulations that the pseudo-phase-space 
profile, the ratio of the spherically averaged DM mass density to the cube of its spheri- 
cally averaged radial velocity dispersion, is well described by a single power-law, Q(r) = 
p(r)/cr r (rf ~ r~ lM , even though neither the density nor the velocity dispersion profiles by 
themselves are [32, 38-40 1. The power law slope is remarkably close to analytic predic- 
tions based on spherical secondary-infall similarity solution [41 1 and their generalization 
Il42l in the inner, virialized regions of halos [43]. Departures from a pure power-law oc- 
cur around the virial radius, close to the location of first shell crossing, where particles 
have not yet fully virialized. Note also that the low velocity dispersion in subhalos leads 
to large fluctuations in local estimates of the phase-space density and thus its spherical 
average does not follow a single power law [32, 43l . 

iii) Substructure 

The numerical resolution achievable by state of the art DM-only simulations has grown 
to the point where it has become possible to follow bound DM structures beyond their 
merging with a larger halo. This has allowed studies of DM substructure, consisting both 
of a population of surviving self -bound subhalos orbiting within the potential of their hosts 
and the debris associated with their tidal stripping and disruption. 

These simulations have been able to probe the subhalo mass function and V max function 
over ~ 5 of magnitude in subhalo mass P4l|45| , and have shown that they are well fit by 
simple power laws, dN/dM ~ M~ L9 and N{> V max ) ~ ^mL- ^ nas even been possible 
to resolve up to 4 levels of the sub-substructure hierarchy [45 1, but the statistics are 
currently not sufficient to quantify sub-substructure scaling laws. 

As with isolated halos, determining subhalo density profiles and concentrations is of 

great interest Il44ll45l . At current resolutions, subhalo density profiles appear to be well fit 
by both NFW or Einasto profiles, and there is some evidence for a radial scaling of subhalo 
concentration, with higher concentrations for subhalos closer to the center [45, 46 1. The 
latter effect is likely due to a combination of stronger tidal forces and the earlier collapse 
times of subhalos found close to the center [47]. The spatial distribution of subhalos 
within their hosts appears to be "anti-biased" with respect to the host's mass distribution 
ll2^IZ7ll431l46ll48ll49l . meaning that the subhalo density normalized by the host's mass 
density profile decreases towards the center. The degree of this anti-bias depends on how 
the subhalo sample is selected: a current mass-selected sample is more affected by tidal 
stripping, which is stronger closer to the host's center, resulting in a pronounced anti-bias 
than a for sample that is selected by properties unaffected by tidal stripping, like the mass 
before accretion iBOllSD . 

Velocity-space substructure, in the form of tidal streams and debris flow, is another topic 
that has received attention E51 15214541 . DM tidal streams have low configuration space 
density, typically only a few per cent of the underlying host halo density. However, since 
they are significantly colder than the host halo, they have a high phase-space density con- 
trast. Another form of velocity-space substructure may be contributed by a so-called dark 
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disk ll55ll , a component of the DM halo that is co-rotating with the stellar disk, which may 
have been deposited by satellites disrupted in the plane of the Galaxy. 

On smallest scales the fine-grained phase-space structure describes the detailed distri- 
bution of DM in configuration and velocity space. Before the onset of nonlinear structure 
formation CDM was almost uniformly distributed with particles lying on a thin three- 
dimensional hypersurface embedded in six-dimensional phase-space. Due to their col- 
lisionless character DM particles then follow the Vlasov-Poisson equations leading to 
stretching and folding of this initial sheet. Therefore, at later times the velocity distribution 
of DM at a given point in configuration space is a superposition of fine-grained streams of 
different velocities. Furthermore, phase-space fold catastrophes lead to caustics, where the 
DM configuration space density can become very large, many order of magnitudes larger 
than the mean halo density ll56l l57ll . It has only very recently become possible to study 
these effects numerically by extending classical N-body schemes Il58ll59l . 

iv) Local DM 

Lastly, numerical simulations provide expectations regarding the DM distribution at the 
solar circle, ~ 8 kpc from the Galactic Center. Simulations indicate that the local density 
of DM is likely to be quite smooth and uniform [52, 60], since strong tidal forces disrupt 
most subhalos close to the center. The flipside of this coin is that the local neighborhood 
should be crossed by many DM tidal streams [52, 53], cumulatively referred to as debris 
flow ED. 

2. 2. Relevance for Dark Matter Detection 

Predictions from cosmological DM simulations affect nearly all DM detection efforts, in a 
variety of ways. In the following sections we highlight some of these inter-dependencies, which 
are also summarized in Table Q] 

2.2.1. Astrophysical Probes 

i) Dwarf galaxies 

The abundance of dwarf satellite galaxies orbiting our Milky Way and M3 1 is potentially 
sensitive to the nature of DM (see §3.2| >. Results from CDM simulations have been used to 
predict how many more ultra-faint dwarf galaxies should be detected, once surveys more 
sensitive than the Sloan Digital Sky Survey (SDSS) and covering the southern hemisphere 
(e.g. DES, Skymapper, Pan-STARRS, LSST) come online [61|. The kinematics of stars 
in the very centers of dwarf spheroidal galaxies have been used to constrain the DM mass 
enclosed within their half-light radius (~ 300 pc) 116211631 and the shape of the DM den- 
sity profile of their host halos, and these measurements have been confronted with the 
predictions from CDM simulations [644 46711 . 

ii) Stellar streams 

The discovery in the SDSS of extended stellar streams 1 68 1, arising from the tidal disrup- 
tion of dwarf galaxies, has provided first-hand evidence for the hierarchical build-up of the 
Milky Way. DM counterparts to these stellar streams are seen in numerical simulations 
Il28l |69l , and raise expectations that many more stellar streams remain to be discovered 
117011 . Additionally, cold stellar streams associated with the disruption of globular clusters 
11711 are promising probes of the dumpiness of the Milky Way's DM halo, since interac- 
tions with passing subhalos should produce kinks and gaps the stream Il72ll73l . 
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Table 1 : A matrix showing which predictions from numerical DM simulations affect which astrophysical probes, indirect, 
and direct DM detection efforts, and vice versa. 



iii) Gravitational lensing 

Gravitational lensing provides important probes of DM on cosmic, cluster, and galactic 
scales that can be compared to the predictions from numerical simulations. We can dis- 
tinguish between weak lensing (see also §2.3 iv), which can be used to tomographically 
map the large scale distribution of DM [74 1 and to measure the total masses and shapes 
of individual halos through cluster and galaxy -galaxy lensing 1751176 1, and strong lensing, 
which can probe the central slope of DM density profiles [77 1 and is sensitive, through 
flux ratio anomalies [78] and potentially time-delay perturbations (79), to the amount of 
DM substructure present in cluster and galaxy halos [80, 81 ]. Recent studies comparing to 
predictions from numerical simulations tend to find that the amount of substructure present 
in DM halos may be insufficient to explain the observed occurence of flux ratio anomalies 
|. However, the effects of intervening line-of-sight structures can be important [ 



2.2.2. Indirect Detection 

Indirect detection of DM refers to the search for the products of pair-annihilations of DM. 
The direct annihilation into two photons is typically loop-suppressed, and so the dominant anni- 
hilation channel is into quarks, gauge (or Higgs) bosons, or directly into leptons. The hadroniza- 
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tion of heavy annihilation products results in gamma ray photons, electrons and positrons, and 
neutrinos. All of these are potentially observable, for example with ground based Atmospheric 
Cerenkov Telescopes (MAGIC, VERITAS, H.E.S.S.) and neutrino detectors (IceCube), balloon- 
borne detectors (ATIC), and space-based satellites (PAMELA, Fermi Gamma-ray Space Tele- 
scope) and experiments (AMS-02 on the International Space Station). In the following we dis- 
cuss some of the possible DM annihilation signatures. 

i) Extra-galactic Diffuse Gamma-ray Background 

The extra-galactic diffuse gamma-ray background (DGRB) refers to the sum-total of all 
gamma-ray radiation produced by DM annihilations throughout cosmic history [86]. The 
amplitude of this signal depends on the large scale distribution of DM in the universe, 
the evolution of the isolated halo mass function, the concentration-mass relationship, and 
of course the density profile. Clumpy substructure may also play an important role, by 
boosting the annihilation luminosity of individual halos (see below). Large scale numeri- 
cal simulations (e.g. Millennium-II |87|) have been used to make predictions for both the 
amplitude of such a signal and the level of its anisotropies [88|, and these have been con- 
fronted with current Fermi data [89 1. The uncertainties of these constraints are dominated 
by the unknown contribution of subhalos below the simulations' resolution limit. 

ii) Galactic Diffuse Gamma-ray Background 

A second DGRB should arise from DM annihilations within the Milky Way's halo, with 
one component stemming from the smooth halo DM and another from clumpy substruc- 
ture. The substructure component is expected to have a shallower angular intensity profile 
than the host halo component [90, 91 1, for two reasons: (i) since it consists of the combined 
emission from a population of subhalos, it should scale with radius as the number density 
of subhalos n su \,(r), rather than as the square of the DM density phost(r) 2 , and (ii) n su \,(r) is 
anti-biased with respect to phostM, with less substructure found close to the host's center. 
DM substructures introduce characteristic anisotropies in the Galactic DGRB |92|, which 
may allow the signal to be disentangled from an astrophysical DGRB. 

The detectability of the Galactic DGRB from DM annihilation thus depends on the abun- 
dance of substructure, its internal structure (concentrations and density profiles), and its 
spatial distribution within the host halo. If the substructure contribution remains sub- 
dominant, the shape of the Milky Way's DM halo may determine the large-scale angular 
variations of the signal. 

iii) Galaxy Clusters 

Galaxy clusters are the most massive gravitationally bound systems in the universe, and 
thus have long been considered good targets for indirect detection searches 11931 . The de- 
tectability of an annihilation signal from clusters relies on a substantial cross section boost 
(either from substructure or Sommerfeld enhancement, see below) |94, 95 1, and the re- 
sulting emission would likely be extended. A difficulty is that any gamma-ray signal from 
annihilation has to compete with the cosmic ray induced gamma-ray emission 11961 . Never- 
theless, H.E.S.S. ||97l and Fermi l95ll data have been able to constrain DM parameters, and 
there is tentative evidence for a ~ 130 GeV line signal from a subset of the most promising 
cluster targets 11981 . 



iv) The Galactic Center 
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The most prominent DM annihilation signal is thought to arise from the Galactic Center 
(GC) [99 1, given its proximity (~ 8 kpc) and the expected high DM density there. Un- 
fortunately, the GC is also an extraordinarily astrophysically active place [ 100 1 hosting 
numerous SN remnants, pulsars, X-ray binaries, and other high-energy sources, not to 
mention a super-massive black hole. Although these astrophysical foregrounds encumber 
DM searches directed towards the GC, it nevertheless has remained a popular target for 
indirect detection efforts. In fact, several gamma-ray "excesses" and anomalies from the 
GC have been reported I101H1041 . including the recent intriguing reports of a gamma-ray 
line at ~ 130 GeV 1 105-108 1. It is too early to confidently claim a detection of DM an- 
nihilation for any of these signals, and additional data will be necessary before statistical 
fluctuations, instrumental effects, or astrophysical sources can be ruled out. 

The strength of the GC DM annihilation signal depends sensitively on the shape of the 
Milky Way host halo's DM density profile at scales that are currently not well resolved in 
numerical simulations. Predictions thus rely on extrapolations of fitting function that have 
been calibrated at larger radii (several hundred of pc) to a small number of high resolution 
simulations OT1 [32l l44l . Furthermore, the gravitational potential is baryon dominated at 
the GC, and one must thus account for modifications of the DM density profile due to 



baryonic physics. As discussed in more detail below (f 3.4 1, these processes may lead to 
either a steepening or a flattening of the profile, and may even displace the location of 
maximum DM density from the dynamical center. 

) Milky Way Dwarf Galaxy Satellites 

The most DM dominated objects known are the Milky Way dwarf spheroidal satellite 
galaxies, which have mass-to-light ratios of up to a 1000 M109I . They are thus natural can- 
didates for indirect detection searches [ 1 10 1 . Since their distances are fairly well known, 
the detectability of their DM annihilation signal is determined by the mass, concentration, 
and density profile of their DM host halos. For many of these systems, stellar kinematic 
data has provided tight constraints on the enclosed mass within the stellar half-light ra- 
dius ll63l . under assumptions of equilibrium and spherical symmetry. Fermi data from 
individual and stacked dwarf galaxies have provided some of the most stringent limits on 
the DM annihilation cross section, extending to below 3 x 10~ 2f> cm 3 s for a DM parti- 
cle mass of ~ 40 GeV and annihilation into pure bb HI 1 II . but these limits assume cuspy 
NFW-like DM density profiles and may be significantly weakened if baryonic processes or 
departures from the CDM assumption result in a flatter profile than inferred from DM-only 
simulations. 

) DarkSubhalos 

The vast majority of subhalos predicted in CDM cosmology are expected to be completely 
dark, since their masses are too low to have allowed gas to cool and form stars 11 121 . 
Individual dark subhalos may be promising sources for indirect detection searches, and 
results from high-resolution simulations have been used to quantify their detectability 19T1 
II 131 II 14l . The Fermi point source catalog contains hundreds of "unassociated" sources 
without identified astrophysical counterparts 111 151 . and it is possible that DM subhalos 
may be lurking among them |116|. Very nearby sources could appear as faint, spatially 
extended gamma-ray sources to Fermi ||9TI . and it may even one day be possible to measure 
proper motions of very nearby subhalos 11 1 1 71 . 
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Once again, these results are highly uncertain due to insufficient knowledge of the abun- 
dance, spatial distribution, and luminosity-mass relation of subhalos on scales below the 
simulations' resolution limit, as well as their ability to survive interactions with the Galac- 
tic disk. 

vii) Local Anti-matter 

DM annihilations in the local neighborhood would produce high energy positrons and 
anti-protons, either through direct annihilation into leptons {e~ e + , /j~ /u + , t~t + ) or via the 
hadronization and decay of other annihilation products. These high energy anti-particles 
might be detectable as an excess over astrophysical cosmic ray backgrounds, and have 
been searched for by the Fermi Oil, H.E.S.S. ESI, PAMELA lIT20ll . ATIC-2 EH, and 
AMS-02 11221 experiments, among others. Owing to energy losses from inverse Compton 
and synchrotron radiation, the propagation distance for positrons is short (~ 1 kpc), and 
thus any injection of positrons from DM annihilations would have to originate from the 
local neighborhood. The expected contribution from DM annihilations hence depends 
on the local density of the Milky Way halo at 8 kpc. The presence of subhalos within 
~ 1 kpc of Earth could affect both the amplitude of this signal and its energy spectrum 
H123II . Improved numerical simulations with higher resolution and accounting forbaryonic 
physics effects will be necessary to better characterize the role of local DM annihilations 
in the high energy cosmic ray spectrum. 

viii) Neutrinos from Earth & Sun 

DM annihilations occurring in the center of the Sun or Earth could produce high energy 
neutrinos that may be observable with neutrino observatories like AMANDA [124| and 
IceCube 11251 . DM particles can be captured by the Sun and Earth through elastic scat- 
tering off of heavy nuclei 111261 . Subsequent scatterings then thermalize the population 
of bound DM particles, and an equilibrium is established between annihilations and cap- 
ture. The strength of the signal depends on the local DM density, but additionally also on 
the speed distribution of incident particles, since lower speed particles are easier to cap- 
ture II 1271 . Rates are thus especially sensitive to the presence of a "dark disk" component 
H128L which can result in a marked increase in the fraction of DM particles traveling at 
low speeds with respect to the solar system. 

ix) Substructure Boost Factors 

Given that the smallest collapsed structures in WIMP CDM are expected to be ~ 10 -12 - 
10~ 3 M Q 11291 [T30I, even the highest resolution numerical simulations can only resolve a 
small fraction of the expected substructure hierarchy. Since annihilation rates scale with 
the square of the density p and (p 2 ) > (p) 2 , any unresolved small-scale dumpiness should 
result in a boost of the annihilation rate. Halo annihilation rates calculated from aver- 
age density profiles (like NFW or Einasto), or even directly from the simulated particle 
distribution, could well underestimate the true luminosity by several orders of magnitude. 
This so-called substructure boost factor has been invoked to motivate effective annihilation 
cross sections orders of magnitude larger than the thermal relic value fe.g. l94l fl 3 1141 34B . 

In the following we discuss two important facts about substructure boost factors that are 
perhaps not as widely appreciated as they should be: 



(a) There is no one single boost factor. 
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Figure 2: An extrapolation of the subhalo contribution to the total luminosity to masses far below the simulation's 
resolution limit. Depending on what one assumes for the concentration-mass relation, one can get very different total 
substructure boost factors. Extrapolations from the high-mass behavior seen in simulations (red dashed) or assuming 
a constant power law concentration-mass relationship (green) are unlikely to hold at masses below ~ 1 M Q (visually 
indicated with thin faint lines). 



The expected substructure boost depends on the distance from the halo center, with 
results from state of the art simulations implying very little (or no) boost at the Galac- 
tic Center, possibly 0(1) in the local neighborhood, and perhaps as large as 100 - 
1000 for the total luminosity of a halo El Ejl ED Q33 EH- As a result a differ- 
ent boost factor applies to spatially extended sources (Galactic DGRB, MW satellite 
galaxies, dark subhalos) than for unresolved sources (distant halos, extra-galactic 
DGRB), and similarly a gamma-ray boost factor may not be the same as those for 
positron or anti -proton production [136]. Furthermore, if a significant fraction of the 
mean density at a given radius is locked up in substructure, then properly accounting 
for the substructure boost will actually lower the smooth density contribution to the 
luminosity 11 141 . further reducing the contrast between the outer regions of a halo 
and its center. The total halo luminosity boost likely depends on the mass of the 
halo, since numerical simulations indicate a roughly equal contribution from every 
decade of substructure mass, and larger mass host halos contain more decades of 
substructure mass ll94l . 

(b) Substructure boosts depend sensitively on subhalo properties many orders of magni- 
tude below the resolution limit of state of the art simulations. 

One approach to estimating the full substructure boost is to stay as close as possible 
to the results from ultra-high-resolution numerical simulations like Via Lactea II and 
Aquarius, by fitting the luminosity boost from all subhalos with mass greater than 
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M m \ n , B(M m i n ) = L(> M m i n )/L smoot h, to a power law of M m ; n over the 4-5 decades of 
substructure mass that are currently resolved in the simulations, and then extrapolat- 
ing this power law down to the free-streaming cutoff scale. This approach was taken, 
for example, by Springel et al. (2008) lfl35l . who found B(M min ) ~ M m °f 26 , and 
inferred a total boost factor for a Milky-Way-like halo of 230 for M min = 10 -6 M Q . 
Another approach is to use the numerical simulation results only to constrain the mass 
function of subhalos, which is measured to be a power law, dnjdM^h ~ M" ub with 
logarithmic slope a ^ -1.9 B4l|45l , and to use an analytical approach to determine 
the subhalo luminosity-mass relation down to the smallest mass halos ||90l |9T1 [T37l . 
The luminosity of a subhalo of mass M is completely determined by its concentra- 
tion c, L/M ~ c 3 //(c), where /(c) depends on the shape of the density profile: for an 
NFW profile, L/M scales approximately as c 2 24 ; for an Einasto profile, as c 2 ' 46 H 1 381 . 
The subhalo annihilation luminosity-mass relation is then completely determined by 
the concentration-mass relation. Again, one may choose to use a simple power law 
relation, for example c(M) ~ M~° • , which well describes the concentration-mass 
relation of Galactic scale halos [33 1. Alternatively one may choose a model in which 
the concentration of a halo reflects the mean density of the universe at its typical col- 
lapse time, as in the analytical model of Bullock et al. (2001) l33l . In this case, the 
concentration-mass relation is not a simple power law, but instead rolls over at low 
masses, and concentrations asymptotically become independent of mass. A com- 
parison of the three approaches discussed so far is shown in Fig. |2j which demon- 
strates how sensitively the total halo boost factor depends on assumptions about the 
small scale behavior of subhalo luminosities. Depending on what one assumes for 
the concentration-mass relation, the total boost of a Milky Way halo ranges from 
3 to 300 (for = 1O~ 6 M ). Note that these three different approaches are not 
all equally likely to apply in reality. Simple extrapolations from the high-mass be- 
havior observed in simulations or assuming a simple power law concentration-mass 
relation are inconsistent with expectation from theoretical models of CDM structure 
formation. Microhalo simulations find concentrations of the smallest and earliest col- 
lapsing DM halos that are incompatible with a single power law c(M) over the full 
substructure hierarchy 11139111411 . Consequently, substructure boost factors much in 
excess of ~ 10 are unlikely to apply in nature. 

A third approach, employed by Kamionkowski et al. (2010) [60|, it to model the vol- 
umetric probability function of density fluctuations, calibrate it to a high resolution 
simulation (Via Lactea II) as a function of halo-centric radius, and then to integrate 
this PDF to obtain an estimate of the boost factor as a function of radius. This ap- 
proach also results in a modest total boost factor for a galaxy-scale halo of O(10). 

Lastly, the annihilation luminosity can in principle also be enhanced by caustics in the 
fine-grained substructure, however recent numerical studies of this caustic boost find only 
percent level increases due to very efficient mixing in phase-space l57l . 

x) Sommerfeld Boost Factor 

A second type of boost factor is of particle physics nature. When the mass of the force 
carrier particle is sufficiently lighter than the DM particle, the so-called Sommerfeld ef- 
fect, a non-perturbative correction for long range attractive forces, can lead to a velocity- 
dependent enhancement in the annihilation cross section 1142111431 . Instead of the usual 
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(crv) = constant behavior, with Sommerfeld enhancement the cross section scales as 
(crv) ~ 1/v, or even 1/v 2 at resonances II 1441 . Although the effect typically saturates 
at small velocities (v/c ~ 10~ 4 - 10~ 5 ) owing to the finite range of the interaction, this 
effect may significantly enhance the annihilation rate in subhalos compared to the smooth 
host halo, given the subhalos' lower velocity dispersions H 145441471 . The details depend on 
the predictions of numerical simulations of the velocity structure in the host and subhalos. 
On the fine-grained level of DM the Sommerfeld effect can have interesting implications. 
Whereas in non-Sommerfeld models the largest annihilation signal is expected to occur 
near caustics due to their high density, this situation changes if Sommerfeld enhancement 
processes are invoked. In that case, cold low-velocity dispersion phase-space structures 
are enhanced compared to hotter regions. Liouville's theorem dictates that DM is very hot 
in caustics to preserve the fine-grained phase-space density. Depending on the details of 
the Sommerfeld model this can make fine-grained streams more prominent for annihilation 
radiation than caustics, because streams are very cold due to their potentially low density. 
This can cause the annihilation rate in streams to dominate over the rate of the smooth 
mean density contribution in halos [ 148 1. 

2.2.3. Direct Detection 

Direct detection refers to efforts to detect nuclear recoil signatures produced in rare DM- 
nucleus scattering events in shielded underground detectors. A large number of experiments 
are currently pursuing this goal, and are utilizing a variety of different technologies and target 
materials (see [149| for a review). The expected event rate and recoil spectrum depends on the 
mass of the target nuclei and of the DM particle, on the nature of the interaction (spin-dependent 
vs. spin-independent), the nuclear form factor, the local DM density po, and the Earth frame 
velocity distribution f(v) of incident DM particles. Until recently, most event rate and parameter 
exclusion calculations assumed a simplified model of the local DM distribution, taking the local 
DM density to be 0.3 - 0.4 GeV cm 3 and a Maxwellian speed distribution with a 1-D velocity 
dispersion of cr — 155 km/s (such that the most likely speed vo is equal to the Galactic rotation 
speed at the solar circle, vo = 220 km/s), and an escape speed of 550 km/s. 

In recent years these assumptions have been directly confronted with the predictions from 
high resolution simulations like Via Lactea II and Aquarius. The large number of self-bound 
subhalos found to be orbiting in the Milky Way's potential raises the question of whether one 
might expect large fluctuations in the DM density at the solar radius. If the Earth happened to 
be passing through a subhalo, for example, the local density of DM might significantly exceed 
the mean value at 8 kpc. Analytical calculations [ 150] and direct "measurements" in simulations 
Il52ll60ll indicate that the volumetric probability distribution of over-densities 6 = p/(p) consists 
of a narrow log-normal reflecting variations in the smooth halo density and a !FV((5) ~ S~ 2 power 
law tail extending over several orders of magnitude before steepening to cT 4 at an overdensity 
corresponding to the mean density of the universe at the collapse time of the smallest halos. 
Barring dramatic changes in the abundance and internal properties of subhalos below the simu- 
lations' resolution limit, the normalization of the power law tail at 8 kpc appears to be too low 
to lead to a non-negligible chance of the Earth lying in a substantial overdensity. It thus seems 
safe to use the mean value of the DM density at 8 kpc in direct detection calculation. However, 
what that value is remains uncertain at least at the factor of two level, with recent studies finding 
values ranging from 0.2 to 0.85 GeV cm" 3 lT5TlfT54"l . 

The speed distribution is another matter. DM-only simulations have definitively demon- 
strated that /(v) shows clear departures from a pure Maxwellian ll27l l52l l53l 11551 11561 . with 
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the typical shape instead exhibiting a deficit near the peak and an excess at lower and higher 
velocities. This is a consequence of the non-Gaussian nature of the three velocity components 
aligned with the density ellipsoid, with the major axis component being platykurtic (broader than 
Gaussian) and the other two leptokurtic (narrower) ll52ll . In addition to these coarse departures 
from a Maxwellian, additional small scale structures are often visible in the high speed tail, in 
the form of broad bumps that are stable in both space and time ll52l l54l [T56 1 and occasionally as 
narrow spikes at discrete speeds indicating the presence of a tidal stream or subhalo in the sample 
volume Il52ll53ll . The presence of a strong "dark disk" [55 1 would change the relative proportion 
of high speed and low speed particles, which could affect scattering rates and the recoil spectrum. 

XCO 
f(v)/vdv, where v m i„ is the minimum speed that can 

result in a recoil of energy which for elastic scattering is given by v m i n = -^m^E^flpp- . 
is the mass of the target nucleus and yu = m^mi^jim^ + wdm) is the reduced mass. The smaller 
mpf, the lower the speeds that are required to produce a given recoil energy. This implies that 
experiments with different target materials and different recoil energy sensitivities probe different 
parts of the speed distribution. Likewise, the mass of the DM particle can affect what range of 
speeds an experiment is sensitive to. For very massive particles (hjdm » m^), the experiment 
becomes insensitive to »jdm, but for so-called "vanilla" WIMP DM with moM ~ 100 GeV, the 
current experiments' Eft-thresholds of ~ 10 keV correspond to v m i„ m 150 km/s. The scattering 
rate will thus be dominated by fairly low speed particles near the peak of /(v) and below. In this 
case, experiments will not be able to see effects from the interesting velocity substructures (the 
bumps and occasional spikes) that occur primarily at high speeds. A strong dark disk, on the 
other hand, may boost event rates. 

On the other hand, with inelastic DM 11571 . for which the relation between Er and v m ; n 
can become inverted, or light DM 11581 (otdm ^ 10 GeV) experiments would be sensitive 
only to much higher speed particles, v m i„ > 400 - 500 km/s. In this case, the departures from 
Maxwellian, both global and local ones, would perhaps be more important. They could alter the 
shape and extent of current parameter exclusion curves [53], potentially reconcile some (but not 
all) of the currently discrepant results from different experiments 111591 11601 . change the phase 
and amplitude of the annual modulation signal ll53l and shift it to higher recoil energies 1541 . 
Directionally sensitive experiments 11611 should be especially sensitive to velocity substructure, 
since they typically have high recoil energy thresholds (~ 50 keV), implying a large v m ; n . The 
majority of high recoil energy particles may in fact be coming from a hotspot significantly offset 
from the direction anti-parallel to Earth's motion through the halo [53 1, and debris flow particles 
can result in ring-like features in the arrival direction [54, 162|. 

Besides generic WIMPs, axions provide another interesting DM candidate. Axions were 
introduced to explain the absence in Nature of strong-CP (Charge conjugation and Parity) vio- 
lations 1163L They are expected to be extremely weakly interacting with ordinary particles, so 
that they never were in thermal equilibrium in the early Universe. This implies that axions can 
be very light (jueV range) and nevertheless form a cold (non-relativistic) component of cosmic 
matter. The cosmological axion population formed out of equilibrium as a zero momentum Bose 
condensate leading to a very small velocity dispersion. In the absence of clustering their present 
day velocity dispersion would be negligible (5v ~ 10~ 17 c compared to Sv ~ 10 _10 c for generic 
WIMPs) making them a good CDM candidate. Axions can be detected through their conver- 
sion to microwave photons in the presence of a strong magnetic field 11641 . Galactic axions 
have non-relativistic velocities (J3 = v/c ~ 10~ 3 ) and the axion-to-photon conversion process 
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conserves energy, so that the frequency of converted photons can be written as: 



v.^ + Av. = 241. 8 ( I ^)(l + i^) M Hz (1) 

where m a is the axion mass that lies between 1CT 6 eV/c 2 and 10~ 3 eV/c 2 . 5/zeV axions would 
therefore convert into = 1200 MHz photons with an upward spread of A ~= 2 kHz due to their 
kinetic energy. An advantage of axion detection compared to WIMP searches is the fact that it is 
directly sensitive to the energy rather than to the integral over the velocity distribution. Narrow 
velocity streams can therefore be more easily detected and lead to spikes in the axion detection 
spectrum. In that case the fine-grained structure can be relevant for detection experiments. A low 
number of fine-grained streams could potentially leave a distinct imprint in velocity-sensitive 
detector signals. For example, recent simulations [57 1 predict that the most massive fine-grained 
streams should be observable with axion detectors like ADMX. 

2. 3. Relevance for Dark Energy Studies 

One of the simplest astrophysical observations, galaxy imaging and the measurement of their 
redshifts and angular positions on the sky, has emerged as a very powerful method to explore the 
nature of DE. For instance, from this data baryon acoustic oscillations, redshift-space distortions, 
abundance of galaxy clusters and weak gravitational lensing can be measured, all of which can 
put constrains on the properties of DE. Consequently, several galaxy surveys (e.g. DES, BOSS, 
BigBOSS, LSST, JPAS, Euclid) are planned or underway to exploit this fact. DM numerical 
simulations have been crucial in this process. 

There are three main areas in which DM-only simulations are essential for the cosmological 
exploration of DE. Firstly, DM simulations have been used to quantify and understand the effects 
of various DE models on structure formation in the Universe fe.g. I165H167I . This allows one to 
identify possible ways of constraining or ruling out some DE models. Secondly, DM simulations 
are the most reliable way to assess potential systematic errors on modeling and cosmological 
parameter extraction from different experiments [e.g |168L Since most of the DE probes involve 
complex and nonlinear processes, an accurate modeling of the signal and related uncertainties in a 
given experiment is of paramount importance in the discovery of new physics. Thirdly, numerical 
simulations can be used to construct mock galaxy and cluster catalogs, which help in the design 
and correct interpretation of surveys aiming to constrain the properties of DE [e.g. 169 1 . Since 
current and future surveys cover large solid angles and probe redshifts out to z ~ 2, simulations 
are forced to very large box sizes and high particle counts, in order to model volumes comparable 
to the surveys while simultaneously resolving structure down to the scale of individual galaxies. 

In the following we briefly describe four of the most established probes aiming to constrain 
the properties of DE and highlight the role of numerical simulations in their development. 

i) Baryon Acoustic Oscillations 

Before recombination, the coupling between free electrons and photons via Thomson scat- 
tering resulted in a distinct pattern of oscillations in the baryon and temperature power 
spectra. These BAOs were also imprinted in the late-time total mass field (albeit with 
smaller amplitude) due to gravitational interactions between baryons and DM. Thus BAOs 
are expected to be present in the galaxy power spectrum, and could be used as a cosmic 
standard ruler IfTTI 1 170441721 . Currently, the feature has indeed been detected at high sig- 
nificance in different galaxy surveys: 2dFGRS fP751 . SDSS lfT74ll . WiggleZ ifTTBI . 6dfGS 
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HI 761 and BOSS M177L placing constrains on cosmological parameters [e.g. |175l 11781 . 
Future surveys are expected to significantly tighten these constraints, in particular those 
on the DE equation of state, which could rule out some DE candidates. 

Large-scale DM simulations played an important role in this development. They showed 
that BAOs survive the diffusing effects of nonlinear evolution and of galaxy peculiar ve- 
locities, and that they should be detectable in biased tracers of the density field I179H1821 . 
At the same time, numerical simulations unveiled significant distortions in the shape of the 
acoustic oscillations due to these effects [180-185 1 which would lead to a systematic error 
on measurements of the acoustic scale. However, recently methods to remove this bias 
have been proposed [186-188 ], and numerical simulations have been used to test their va- 
lidity and performance. In the future, specially-designed numerical simulations will help 
us to understand and model better the impact of structure and galaxy formation on the 
observed BAO signal in the clustering of galaxies. 

) Redshift space distortions 

The redshift of a galaxy not only contains information about its cosmological distance, but 
also about its peculiar velocity. This difference between angular positions and redshifts 
creates an anisotropy in the observed two-dimensional galaxy correlation function that 
can be used to establish the relation between density and velocity fields in the Universe 
[ 1 89 ] . These redshift-space distortions (RSD) have been historically used to constrain the 
value of the matter density of the Universe 1 1901 11911 , but recently have also been em- 
ployed to constrain the gravity law 111751 119214194L Current measurements are consistent 
with General Relativity (GR), but future surveys are expected to significantly improve the 
constraints. Hypothetical departures from GR could be related to the DE model. 

However, extracting this information is not trivial. Numerical simulations have shown that 
linear perturbation applies only on extremely large scales (> 50 Mpc/h) 1 1831 11951 11961 . 
Quasi-linear corrections and the so-called "finger-of-God" (FoG) cannot be neglected on 
smaller scales. In particular, FoGs are produced by the motions of galaxies inside DM 
halos (whose velocity is comparable to bulk motions produced by large-scale density per- 
turbations), which introduces a strong damping in the galaxy clustering along the line of 
sight. 

Given the increasing number of Fourier modes, it is desirable to use RSD down to the 
smallest possible scales. Observations are usually modeled as a combination of linear 
theory expectations plus a damping to account for FoG [e.g. |192| . However, numerical 
simulations have highlighted the pitfalls of this approach and the systematic error that it 
would introduce for future surveys [ 1681 119311197 198 1. This finding has fueled the devel- 
opment of more accurate new estimators which can robustly use the clustering information 
at smaller scales [ 199 200 1 . This is another example of the importance of DM numerical 
simulations in the optimal exploitation of observational datasets. 

) Abundance of galaxy clusters 

The position of galaxies in an optical survey can also be used to identify galaxy clusters. 
The number of these objects a function of their mass is of great interest because it contains 
information about the underlying probability distribution function of density perturbations 
in the Universe. The evolution of the mass function on group and cluster scales has in- 
deed been used to measure cosmological parameters, helping to break degeneracies in the 
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constraints from other probes [201-206]. In addition, the highest mass end is sensitive 
to primordial non-Gaussianities and early DE, and thus the detection of massive, high- 
redshift clusters has been used as evidence of their existence [e.g. 207-213 1. However, 
uncertainties in the mass estimation of clusters and the respective Eddington bias have 
seriously hampered these measurements B2141 12151 . 

The most relevant aspect of DM simulations for this cosmological probe is the halo mass 
function and its dependence on cosmology [e.g. 23 1. This prediction is usually parametrized 
in terms of the linearly extrapolated variance of the underlying DM field 120-221 12161 . 
However, recently numerical simulations have shown evidence for dependencies with 
other parameters 12171 . This finding could seriously limit the maximum performance of 
current approaches to cosmological parameter extraction using clusters. For this reason, 
in the future numerical simulations will probably play a direct role in the modeling of the 
observed abundance of clusters. 

Another important aspect in this probe is in the characterization of the performance of 
cluster finder algorithms. DM-only simulations can quantifying the impact of projection 
effects, misidentification of the cluster's center, and false detections [218 -220], as well as 
the relation between cluster mass and observed richness or weak lensing signal. Hydro- 
dynamical simulations have an analogous role for experiments using X-rays or thermal 
Sunyaev-Zeldovich detected clusters. 

iv) Weak lensing 

The light from high-redshift galaxies gets distorted by intervening mass before it reaches 
us. Deep gravitational potentials cause large distortions which can split the image of a 
galaxy into multiple lensed images, in a phenomenon known as strong gravitational lens- 
ing. This has been discussed in 2.2.1 and can be used to probe the mass of DM halos and 
even the law of gravity and the nature of DM. Smaller distortions in the shape (and size) of 
background galaxies caused by the cosmic web are known as weak gravitational lensing. 
These changes in the properties of galaxies can be related to integrals of the DM mass 
power spectrum and thus they can be used to reconstruct the full three-dimensional DM 
density field. This allows measurements of the growth of structure, which can be used to 
constrain DE and modified gravity 12211 . 

The shear of galaxies has been detected statistically in many surveys [222 -228 ] and has 
been used to place constraints on cosmological parameters [229 -231 ]. For these, the main 
ingredient is the dependence of the nonlinear DM correlation function on cosmology, 
which is normally taken from fitting formulae calibrated using predictions of DM sim- 
ulations [232 -236]. However, weak lensing measurements are affected by many sources 
of systematic errors: most importantly the intrinsic alignments in the shape of galaxies 
caused by tidal forces 12371 , as well as the PSF ellipticity caused by atmospheric dis- 
tortions M238I . among others. Over the last few years extensive studies of these effects 
have been carried out, with DM simulations helping to create synthetic data as well as 
constraining the impact and magnitude of intrinsic alignments. 

The next generation of surveys are expected to be able to reduce systematic effects drasti- 
cally, and thus will require high-precision predictions of the nonlinear DM power spectrum 
down to small scales. Since perturbation theory approaches can provide a reasonable de- 
scription only in the mildly nonlinear regime, the necessary predictions and modeling of 
data will have to rely on DM numerical simulations, either directly or via emulators. 
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3. Current State of the Art 



In this section we present a late 2012 snapshot of the state of the art of cosmological numer- 
ical simulations, with a focus on runs with particular relevance to the DM and DE problems. We 



first discuss DM-only simulations (5 3.1 1, which are mature, mostly computational resource lim- 



ited, and have been pushed to extremely high resolution, and then DM+hydro simulations (§ 3.4 1, 
which are algorithmic ally more challenging, less well developed, limited to lower resolution, and 
do not yet produce robust or even converged results. 



3.1. Dissipationless Dark-Matter-only Simulations 
3.1.1. Numerical Techniques and Codes 

Pure DM simulations take the ansatz of completely neglecting any dissipational baryonic 
physics and treat all matter as collisionless DM. The density field is sampled with discrete "N- 
body" particles, whose gravitational evolution is governed by the Poisson-Vlasov equations in 
a coordinate system that is co-moving with the mean expansion of the universe. The effects of 
DE are generally confined to the expansion history, i.e. the translation between cosmic time and 
expansion scale factor [239-241]. Many different techniques have been developed to solve this 
set of equations, and we refer the reader to [ 242-246 1 for extensive discussion. For the present 
purpose, it suffices to briefly describe two of the major techniques in use today. 

One is the so-called tree code 112471 . in which the particle distribution is organized in a hi- 
erarchical tree structure. Contributions to the gravitational potential from distant particles are 
approximated by the lowest order terms in a multipole expansion of the mass distribution at a 
coarse level of the tree. If accuracy requirements demand it, the tree is "opened" to a higher level 
and a more detailed particle distribution is accounted for. This method reduces the computational 
complexity of the N-body problem from 0(N 2 ) to 0(N log N), with a well controlled error. A 
further improvement to 0{N) scaling is possible through the use of the Fast Multipole Method 
(FFM) [248 1, in which forces are calculated between two tree nodes rather than between individ- 
ual particles and nodes. In order to avoid unphysical hard scatterings between nearby particles 
(which are just tracers of an underlying continuous density field), gravitational interactions are 
"softened" on small scales |249|, typically either with a Plummer or a cubic spline kernel. The 
force resolution of this method is then given by the softening length e so f t , which in DM-only 
simulations is usually kept constant in co-moving coordinates. Pkdgrav2 |244| is a prominent 
example of a pure tree code, and it uses FFM. 

The second commonly used N-body technique is the adaptive particle-mesh (PM) method, 
in which the particles are deposited onto an regular mesh to produce a density field. The mesh 
structure is often adaptively refined in high density regions demanding increased force accuracy. 
The gravitational potential is obtained via Fourier transform on the periodic root grid (coarsest 
level), and a multi-grid relaxation technique is used to evaluated it on the refined grids. This 
method also achieves 0(N log N) scaling, but here N refers to the number of mesh cells, which 
is typically taken to be 2 3 times the number of particles. No explicit force softening is neces- 
sary, since particles interact with each other through a mean field not individually, and the force 
resolution is effectively given by the cell spacing of the most refined mesh. Examples of pure 
adaptive-PM codes are Art B242I and Ramses |250|. 

One of the most widely used cosmological simulation codes is the hybrid Tree-PM code 
Gadget 1245 L which uses the PM method to evaluate long range forces and the tree method for 
short range interactions. The Gotpm code 125 11 is another example of such a hybrid. 
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The choice of gravitational softening length in cosmological simulations with tree codes 
is a contentious issue that has not been truly resolved. The difficulty arises because there are 
conflicting demands on the softening 12521 12531 : on the one hand it is desirable to choose as 
small of a softening as allowed by computational resources (the smaller the softening, the shorter 
the time steps, and the more expensive the simulation), since it represents a distortion of the true 
gravitational potential and leads to overmerging [254]. On the other hand, smaller softening 
results in stronger unphysical two-body relaxation effects, which can cause spurious heating as 
well as artificial fragmentation [255, 256]. Some studies have advocated for softening lengths no 
smaller than the mean particle separation in the initial conditions [255, 256 1, while others have 
argued that it is sufficient to choose a softening that suppresses unphysical discreteness effects in 
collapsed region 112571 12581 . Cosmological zoom-in simulations (see 5 3.1.21 generally employ 
softenings significantly below the mean initial condition separation of high-resolution particles 
(e.g. ranging from 0.006 to 0.02 for the zoom-ins in Table [2J. As an aside, the need to avoid 
two-body relaxation effects is responsible for the slow * N 1 ^ rate of convergence of halo density 
profiles [258 1 . Note that due to the hierarchical nature of collapse in CDM, no matter how many 
particles are used in a simulation, the first structures to collapse are always the smallest halos that 
are resolved with only a small number of particles and hence susceptible to 2-body relaxation 
effects. 



3.1.2. Simulations 

Cosmological DM-only simulations can be divided into two types: (i) full-box and (ii) 
zoom-in simulations. The former resolve the entire computational domain with a single par- 
ticle mass and force resolution, and typically cover a representative volume of the universe, with 
box sizes ranging from ~ 100 Mpc to tens of Gpc. They are generally focused on resolving the 
large scale structure of the universe and are most useful for statistical studies of DM halos. We 
refer to this class as cosmic scale simulations. 

In the zoom-in class, simulations forgo capturing a representative fraction of the universe, 
and instead focus all available computational resources on one halo of interest, resolving its 
internal structure and substructure at the highest possible resolution. In order to achieve this 
goal, these simulations make use of nested initial conditions, in which the great majority of the 
computational domain is sampled only with very coarse resolution (large particle masses and 
force softenings), but a small volume containing an object of interest is resolved with much 
higher resolution. We distinguish between cluster scale simulations, in which the object of 
interest is a single galaxy cluster (10 14 - > 10 15 M Q ), and galactic scale simulations, which 
zoom in on a single galactic halo (< few x 10 12 M G ). In both cases the halo of interest is typically 
identified in the z = output of a lower resolution full-box simulation. The particles contained 
within some multiple (usually 3-5) of its virial radius are traced back to the initial conditions, and 
the low resolution particles in this Lagrangian volume are replaced with a nested set of higher 
resolution (lower mass) particles. The phases and amplitudes of the large wavelength Fourier 
modes used to calculate the initial condition displacements of these higher resolution particles are 
kept the same as in the coarse realization, but additional random modes are introduced at smaller 
wavelengths. This process ensures that the large scale matter distribution remains identical to the 
coarse simulation, but with greatly enhanced power in small scale substructure fluctuations. The 
technical details of this procedure are discussed in detail in the literature 1264- 42671 . 

In the following we briefly go over each of these three classes of simulations (cosmic, cluster, 
and galactic), highlighting both state of the art achievements as well as limitations and shortcom- 
ings. Over the last decade progress in this field has been driven by advances in computing 
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For AMR simulations (Ramses, Art) e so f t refers to the highest resolution cell width. 
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Aquarius A- 1 

GHalo 
Via Lactea II 


Gadget-3 
Pkdgrav2 
Pkdgrav2 


5.9 
3.89 
4.86 


4.3 x 10 9 
2.1 x 10 9 
1.0 x 10 9 


1.7 x 10 3 

1.0 x 10 3 

4.1 x 10 3 


20.5 
61.0 
40.0 


82 
43 
13 


El 
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Table 2: Current state of the art in DM-only simulations on cosmic, cluster, and galactic scale, ordered by number of 
simulation particles. Lhires is a proxy for the size of the high-resolution region in zoom-in simulations, and is defined to 
be equal to the size of a cube at mean density enclosing all high resolution particles. N^^ ub is the number of halos 
in the box (Cosmic) or subhalos within (Cluster and Galactic) with at least 100 particles at z = 0. In some cases 
(DEUS FUR, Horizon-411) mass functions have not been published, and so we estimated N^™ 13 from a Sheth & Tormen 
1 19 1 mass function fit. 

technology and available resources at national supercomputing facilities, and were aided by the 
algorithmic developments discussed in the previous section. The simulations discussed below 
all required multiple millions of CPU-hours on thousands of processors, and required terabytes 
of memory and petabytes of disk storage. Some of the characteristics of these simulations are 
summarized in Tables |2]and[3] and visualizations for a subset are shown in Fig. [3] 

i) Cosmic scale 

In this class the state of the art has reached > 10 billion particle simulations, with the current 
record holder (in terms of particle number), the recently completed DEUS Full Universe 
Simulation (FUR) 1125911 , utilizing more than half a trillion particles in a 21 hr l Gpc box, 
which corresponds to the entire observable universe. It was run with a modified version 
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12 
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3 


15 
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Germany 


3.5 


1024 


3 


45 
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1 


60 
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Cray XT4 


Oak Ridge National Lab 


USA 
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20 



Table 3: Supercomputers and computational resources utilized for each simulation. 



of the Ramses code and took about 10 million CPU-hours using 38,016 MPI tasks on the 
European supercomputer Curie. With a particle mass of 1.2 x 10 12 Mq and force resolution 
of 40h~'kpc, DEUS FUR cannot resolve individual galaxies and thus is limited to studying 
the large scale distribution of matter in the universe. The main driver of this simulation has 
been to quantify the imprint that DE leaves on cosmic structures (e.g. BAO), and how the 
nature of DE may be inferred from observations of large scale structure. In addition to a 
standard ACDM (w = -1) run, two more FUR simulations at the same resolution but with 
different DE models (w = -0.87 Ratra-Peebles quintessence and w = -1.2 phantom fluid) 
have recently been completed. 

Currently the only calculation able to simultaneously resolve scales relevant for BAO detec- 
tion as well as those DM halos and subhalos expected to host galaxies to be seen in future 
surveys is the Millennium-XXL simulation. This simulation uses slightly fewer particles 
(303 billion) than DEUS FUR to represent the mass field in the Universe, but it has almost 
200 times better mass resolution due to its smaller computational domain. It was run during 
the summer of 2010 at the Julich Supercomputer Centre in Germany using 12,288 CPUs 
using a memory-efficient version of the Gadget-3 code. The main goal of this simulation is 
to explore the impact of galaxy formation physics on cosmological probes, in particular for 
BAO detection and redshift-space distortion tests. 

On considerably smaller but still cosmic scales, two of the most prominent simulations 
are the Millennium-II and the Bolshoi simulations. Millennium-II, a Gadget-3 simulation, 
has 10 billion particles in a lOOfr'Mpc box, for a particle mass of 6.9 x 10 6 M o . It cost 
1 .4 million CPU-hours on an IBM Power-6 supercomputer at the Max-Planck Computing 
Center in Garching, Germany. Bolshoi, an Art simulation, uses 8.6 billion particles in a 
250h~'Mpc box, giving a particle mass of 1.4 x 10 8 M , and required 6 million CPU-hours 
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on the Pleiades supercomputer at NASA Ames. Both simulations have a force resolution 
of ltr'kpc. Although Bolshoi has 20 times poorer mass resolution, it covers 16 times 
more volume than Millennium-II. One additional difference between the two is the choice 
of cosmological parameters, with Millennium-II employing values inspired by the first year 
WMAP results (O m = 0.25, Q. A = 0.75, h = 0.73, cr 8 = 0.9, and n s = 1), which for cr 8 
and n s are more than 3cr discrepant with the more recent WMAP 5-year and 7-year results, 
while Bolshoi used values (Q,„ = 0.27, Q A = 0.73, h = 0.70, cr 8 = 0.82, and n s = 0.95) 
that are consistent with the more recent measurements^ For both cases, the mass and force 
resolution is sufficient to resolve some of the internal (sub-)structure of Milky Way-like 
halos, while at the same time capturing a large enough sample of such galaxies (~ 5000 in 
Millennium-II, ~ 90,000 in Bolshoi) to enable statistical studies. These simulations have 
provided precise and robust results on DM halo statistics like the mass function, subhalo 
abundance, mass and environment dependence of collapse times, and spatial correlation 
functions, spanning a wide range of scales, from dwarf galaxy halos to rich galaxy clusters. 

ii) Cluster scale 

About ten years ago, cluster scale DM-only simulation were leading the effort to study the 
properties of individual DM halos and the abundance and properties of their substructure 
11271 12701. With around ten million high resolution particles, these simulations resolved 
thousands of subhalos and established important substructure scaling relations. DM sub- 
structure studies then shifted focus to the Galactic scale (see below), and until very recently 
cluster simulations had not pushed into the billion particle regime. The Phoenix simulation 
suite {263 1 has now changed that, with their highest resolution Gadget-3 simulation em- 
ploying 4.1 billion particles to resolve a 6.6 x 10 14 h~ l M cluster, and identifying a total of 
almost 200,000 individual subhalos (60,000 with more than 100 particles). This simulation 
was run on 1024 cores of the DeepComp 7000 supercomputer of the Chinese Academy of 
Science and cost 1.9 million CPU-hours. Additional simulations of rare and dynamically 
young objects like galaxy clusters will help to clarify to what degree the internal structure 
of DM halos and substructure scaling laws are universal and self-similar. 

iii) Galactic scale 

On Galactic scales the three flag-ship simulations are Via Lactea II PHI . Aquarius A-l |45l . 
and GHalo (i32J, in chronological order. Via Lactea II (1.5 million core-hours on Oak Ridge 
National Lab's Jaguar), was the first simulation to use over a billion high resolution particles 
to resolve a single halo, Aquarius A-l (3.5 million core-hours at the Leibniz Rechenzentrum 
in Garching, Germany) the first to have over a billion particles within the virial volume of 
the halo, and GHalo (2 million core-hours on Marenostrum at the Barcelona Supercomput- 
ing Center, Spain) is currently the simulation with the highest mass resolution. With particle 
masses ranging from 1000 to 4100 M Q and force resolutions from 20 to 60 pc, these sim- 
ulations are able to resolve in unprecedented detail the formation and accretion history of 
Milky Way-sized DM halos (M « 10 I2 M o ), their inner density profiles, and the properties 
and survival of stripped subhalo cores, as well as tidal debris orbiting within these systems. 
Density profiles have been reliably measured to ~ 100 pc, and the substructure hierarchy 
is resolved over five decades in mass, down to ~ 10 5 M Q subhalos. The Aquarius project 



3 Results from the Millennium simulations have been rescaled to the latest set of cosmological parameters 1 268 269 1. 
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Figure 3: Visualizations of state of the art simulations on cosmic (Millennium-XXL [220 ], upper left), cluster (Phoenix 
A-l 12631 . upper right), and galactic scale (Aquarius A-l 1 45 J, lower left, and GHalo 1321 , lower right). 



simulated an additional five halos at somewhat lower resolution (particle mass w 10 M Q ), 
which has enabled a valuable initial assessment of halo-to-halo scatter. 

Even though the simulations were run with different codes (Gadget-3 for Aquarius, Pkd- 
grav2 for Via Lactea II and GHalo) and used somewhat different cosmological parameters 
(most notably (<x 8 , n s ) - (0.9, 1 .0) and (0.74, 0.95), respectively), the numerical results agree 
remarkably well with each other when scaled by the mass of the simulated host halo. Some 
disagreements persist, however, in the interpretation of these results, for example in the as- 
sessment of the relative detectability of the Galactic DGRB indirect detection signal and 
that from individual subhalos ||9T1|135L and in the self-similarity of the (sub-)substructure 
population Il44ll45ll . 
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As will have become clear from the previous sections, these state of the art DM-only sim- 
ulations on cosmic, cluster, and galactic scales require truly massive computational efforts (see 
Table [3]). Note that computational demands do not scale solely with N-body particle count, but 
also sensitively depend on the degree to which the simulations resolve small scale structure and 
non-linear clustering, mostly because more time steps are required. Single halo zoom-in sim- 
ulations also require more communication and it is more difficult to balance their memory and 
CPU requirements than for full-box single resolution runs. For this reason a galactic-scale sim- 
ulation like Via Lactea II required about the same number of CPU-hours (about one million) as 
the cosmic-scale Millennium-II run, even though the latter employed 10 times more particles. 
In addition to high computational demands at run time, simulations at this level present enor- 
mous challenges for data transfer, storage, and analysis (see 54.3.1 ii). As detailed nicely in the 
DEUS FUR simulation paper [259], analyzing the simulation often requires computing resources 
comparable to running it. 



3.2. Small scale challenges for Cold Dark Matter 

Predictions from CDM simulations of the large scale distribution of DM, post-processed to 
include mock galaxy populations, agree remarkably well with the observed clustering of galaxies 
measured in modern surveys like the Sloan Digital Sky Survey (SDSS) [18|. Yet at smaller 
scales the agreement between CDM predictions and observations is not as good: the number of 
dwarf satellite galaxies observed to be orbiting our Milky Way (and our nearest neighbor galaxy, 
M31) is less than one would naively infer from the predictions of DM-only simulations in a 
CDM cosmology 1127 II 12721 . The severity of this so-called Missing Satellites Problem has been 
reduced in recent years through the discovery of more than ten new ultra-faint dwarf satellites 
in the SDSS 11091 and references therein], raising the possibility that hundreds more remain 
yet to be discovered [61 1. Nevertheless, reconciling the steep slope of the DM subhalo mass 
function with the shallow faint end of the satellite luminosity functions remains a theoretical and 
computational puzzle. 

A second major challenge to CDM is the Cusp/Core Controversy concerning the central 
slope of DM density profiles in low mass galaxies. Two-dimensional stellar and gaseous kine- 
matic measurements in low surface brightness field galaxies I1641I273H275I . as well as chemo- 
dynamical measurements in at least two Milky Way dwarf satellite galaxies ll65l . imply that the 

slopes of the DM density profiles are shallower than the NFW slope of 1 predicted by CDM 

simulations without baryonsj^] 

Lastly and possibly connecting the previous two concerns, it has recently been pointed out 
B2771 12781 that there may be a problem in the abundance of even the most massive Galactic 
subhalos. Dubbed Too Big To Fail, this problem refers to the inference from stellar kinematic 
data that the central densities of the "classical dwarfs" (bright satellites with luminosities greater 
than 1O 5 L ) are too low to be consistent with inhabiting the most massive subhalos predicted 
in the Via Lactea II and Aquarius simulations. The consequence being that either there exists a 
population of massive subhalos orbiting within the Milky Way's virial volume that have remained 
completely dark and devoid of stars, despite the fact that less massive subhalos manifestly were 
able to form galaxies, or that the Via Lactea II and Aquarius halos are somehow not representative 
of our Milky Way. For example, if the mass of the Milky Way were a factor of two less than in 



4 Cored DM distributions have also been inferred for more massive spiral galaxies le.g. 12761 . but given that these 
systems are strongly baryon dominated, this observation is not commonly considered a major challenge for CDM. 
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these simulations (but still within the range allowed from observational constraints), then the 
number of discrepant (too dense) halos may be small enough to not be a major worry 12791 12801 . 
Alternatively, some process not captured in the DM-only simulations could reduce the central 
densities in the most massive subhalos. 

3.3. Simulations with Departures from Collisionless Cold Dark Matter 

The two assumptions that underlie all DM-only simulations described so far are (i) that DM 
is "cold", meaning that the cutoff in the density fluctuation power spectrum occurs on scales 
far below what is resolved in the simulations, and (ii) that it is collisionless, meaning that the 
only dynamically relevant interactions are gravitational, i.e. that any self-scattering effects are 
negligible. Although these assumptions are theoretically well motivated, holding for example 
for most supersymmetric DM models as well as for axions, they are not a priori requirements. A 
number of studies in the literature have investigated whether departures from the assumptions of 
cold and collisionless DM can provide solutions to the small scale challenges to CDM discussed 
in the previous section. 

In the Warm Dark Matter (WDM) scenario the DM particle exhibits some non-negligible 
thermal velocities at high redshifts, instead of being truly cold. In this case, free streaming in the 
early universe will erase small scale density fluctuations, preventing the formation of low mass 
DM halos. For a WDM particle of mass m x and temperature T x , a cutoff in the power spectrum 
then occurs Ifl2l|281|| at a scale of 



Here T x has been expressed in terms of the temperature of the cosmic neutrino background 
T v , and (T x /T v ) in general is model dependent. For thermal relic WDM particles (e.g. the 
gravitino 02821 ) it can be related to the relic DM density via Q. x h 2 = (T x /T v ) 3 (m x /94eV), but 
non-thermally produced particles (e.g. the sterile neutrino [283 1) can have a wide range of tem- 
peratures. Constraints from the Lyman-a forest limit the mass of a WDM particle to be > 2 - 8 
keV M281I I284 285 1, depending on the details of the particle physics. 

The suppression of small scale power may help to explain the puzzling dearth of Milky Way 
satellite galaxies [286 , 287 1. A secondary effect arising from the lack of small scale structure is 
that the collapse times of halos above the free-streaming cutoff are delayed. This results in lower 
concentrations and reduced central densities, which may help to address the Too Big To Fail 
problem M288L Low concentration halos are also more prone to tidal disruption, which further 
reduces the abundance of low mass objects in WDM halos. Lastly, we mention that WDM halos 
are expected to have central density cores, since the WDM particles' non-zero temperature results 
in a finite phase-space density in the early universe B289L which by Liouville's theorem cannot 
grow during the formation of a halo. For realistic models that are consistent with constraints 
from the Lyman-ff forest, however, it can be shown [290, 291] that phase-space density limited 
cores only occur on very small scales, r core /R v i r < 10~ 3 , far below where there is observational 
evidence for a flattening of the DM density profile. WDM models by themselves thus do not 
appear to be capable of solving the cusp/core controversy. 

There are several technical difficulties associated with numerically simulating WDM mod- 
els. One is that the presence of a cutoff in the power spectrum in the initial conditions gives 
rise to the formation of a large number of spurious halos of purely numerical origin 12921 . An- 
other difficulty is that for sufficiently light WDM particles, small box sizes, or early simulation 




(2) 
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starting times the thermal velocities can become comparable to the bulk flows induced by the 
density fluctuations in the simulation's initial conditions. One should then apply thermal stream- 
ing velocities to the N-body particles, ideally by splitting each particle into 2N sub-particles and 
applying equal and opposite velocities randomly drawn from the primordial velocity distribution 
to each N pairs of sub-particles M293I . In practice, however, thermal streaming velocities can 
usually safely be neglected, unless one is simulating WDM with m x < 1 keV (observationally 
ruled out) or using a boxsize smaller than 1 Mpc. 

The hypothesis of essentially collisionless DM has also been contested. This leads to the 
idea of Self-Interacting Dark Matter (SIDM) ll294i]297l . Initially SIDM models with a con- 
stant scattering cross section were quickly abandoned since those that could solve the small-scale 
CDM problems seemed to violate several astrophysical constraints, such as the observed ellip- 
ticity of the mass distribution of galaxy clusters [298] and the survivability of satellite halos 
112991 . But recently it was pointed out that some of these earlier constraints were overstated, and 
small velocity-independent self-interaction cross sections can have sizable effects on halo profiles 
without violating astrophysical constraints [300|. Also simple ad hoc velocity -dependent cross 
sections of the form l/v a were explored [301 ), yielding encouraging results that however lacked 
a proper underlying particle physics model. More recently it was realized that self-interactions 
through a Yukawa potential can resolve the challenges facing velocity-independent SIDM mod- 
els [ 302 1. The velocity dependence of scattering through the massive mediator of this dark force 
(similar to a screened Coulomb scattering in a plasma) could make scattering important for dwarf 
galaxies with low velocity dispersion, but unimportant at the much higher velocities encountered 
in galaxy clusters. Such models have been explored numerically and it has been shown that they 
can help to resolve some of the small-scale CDM problems through the formation of a central 
density core 13031 . Note that there also exist hybrid models (e.g. Atomic DM lfT3l . see Fig. [TJ, 
in which the DM exhibits both self-interactions and suppression in small scale power. 

3.4. Simulations Including Baryons Physics 

As we have seen, the tension on small scales between dwarf galaxy observations and the 
predictions of DM-only simulations might be an indication that the true properties of the DM 
particle differ from the cold and collisionless assumptions of these simulations. Unfortunately, 
however, the effects of modified DM particle physics can be mimicked by the complicated bary- 
onic physics governing the formation of stars and galaxies inside DM halos. For this reason 
the problem of DM is closely coupled to the problem of galaxy formation, which of course is 
a worthy topic of study in its own right. A survey of the current state of the art in numerical 
simulations of galaxy formation is considerably beyond the scope of this review, and so in the 
following we instead provide a limited overview of recent results with particular pertinence to 
the DM and DE problems. 

3.4.1. Numerical Techniques and Codes 

The basic equations that are solved in cosmological hydrodynamics simulations are the Euler 
equations (conservation of mass, momentum, and energy) governing the flow of an ideal gas, 
coupled gravitationally to the DM sector through a source term in the energy equation and the 
Poisson equation. Neglecting viscosity (ideal gas) is a good assumption on cosmological and 
galactic scales, but ignoring the effects of magnetic fields and radiation less so, and accounting 
for magneto- and radiation-hydrodynamic effects in cosmological galaxy formation simulations 
is an active area of research [e.g. 1304443081 . Results from such studies, however, have not yet 
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been brought to bear on the DM and DE problems, and so we focus here on the simpler pure 
hydrodynamic case. 

Unlike for the purely gravitational N-body problem, where even conceptually quite differ- 
ent solvers (e.g. tree and adaptive PM codes, see 5 3.1.1 1 robustly produce similar results, the 
choice of method with which to treat the hydrodynamics can lead to marked differences in the 
results 009443 15l . Numerous detailed discussions of the different approaches and their relative 
advantages and disadvantages exist in the literature [e.g. 131 1113 13113 161 13 171 . and we only briefly 
summarize the essentials here. 

In general one can distinguish between Eulerian and Lagrangian methods, which discretize 
either space (Eulerian) or mass (Lagrangian). In Smoothed Particle Hydrodynamics (SPH) 
B317L the most commonly used Lagrangian approach^ the fluid flow is followed with parti- 
cles, whose equations of motion are derived from a discretized particle Lagrangian [318], which 
ensures excellent conservation of mass, momentum, energy, entropy, and angular momentum. 
Thermodynamic quantities (density, pressure, etc.) are obtained by smoothing over neighbor- 
ing particles with a particle-dependent smoothing length. Advantages of SPH are that it is au- 
tomatically adaptive, delivering higher resolution in collapsing regions, geometrically flexible, 
inherently Galilean invariant (errors don't depend on bulk flows), computationally inexpensive, 
and that it couples easily to an N-body gravity treatment (as for the DM). It is also often sim- 
pler to implement new physics prescriptions. However, SPH is not without its drawbacks. Its 
Lagrangian nature results in less resolution in lower density environments, and its estimates of 
thermodynamic quantities are noisy on the scale of the smoothing kernel. Artificial viscosity 
must be added in order to inject the entropy generated at shocks and to suppress unphysical os- 
cillations in the states immediately surrounding it. This broadens the discontinuity to several 
smoothing lengths and has a tendency to make the method more dissipative. It has low accuracy 
for contact discontinuities, and as a result suppresses some astrophysically relevant fluid insta- 
bilities and mixing. However, not all of these disadvantages are inherent to the SPH method, and 
several recently proposed modifications have resulted in very promising improvements 113194 - 
322 1 . SPH simulations commonly employ a gravitational force softening length that is fixed in 
physical coordinates below some redshift (see B323I1 ). which results in poorer spatial resolution 
at early times compared to a constant co-moving softening scale and may suppress early star 
formation. The most widely used SPH codes in the galaxy formation field are Gadget 112451 and 
Gasoline 13241 . 

With Adaptive Mesh Refinement (AMR), the most widely used Eulerian hydrodynamics 
approach, the fluid flow is instead discretized in space. Euler's equations are solved on a regular 
mesh, which is adaptively refined in regions requiring higher accuracy 113251 . In finite volume 
methods, conserved quantities (mass, momentum, energy) are stored on the cells of the mesh, and 
their values are updated in a conservative fashion by solving for the fluxes across cell interfaces 
(for a primer, see 113161 ). In the widely used Godunov schemes, fluxes are calculated by con- 
sidering the independent variables to be piecewise constant across each cell, with discontinuities 
at the cell boundaries, and solving the resulting Riemann problem for the characteristic waves 
(shocks, rarefactions, and contact discontinuities) traveling into the neighboring cells 13261 . In 
practice higher order accurate methods (piecewise linear or piecewise parabolic) are commonly 
employed [327-329]. Advantages of AMR are its low noise and high accuracy for shocks, con- 
tact discontinuities, and shear waves, allowing it to capture fluid instabilities with high fidelity, 



Really it is "pseudo-Lagrangian", since shearing flows with distinct internal properties are not followed in a truly 
Lagrangian way on scales below the smoothing kernel |314|. 
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Figure 4: Visualizations of three recent cosmological hydrodynamical galaxy formation simulations. Top row: gas 
surface density at z = 4 in three galaxies (out of ~100 in the box) simulated with the AMR code Enzo 1 336 1. Middle row: 
a series of zooms onto the density field surrounding a z = 2 galaxy, simulated with the moving mesh code Arepo 13141 . 
Bottom: Optical and UV composite images of Eris, a Milky Way-like galaxy simulated with the SPH code Gasoline 
(3371 . 



and its full control over where to place high resolution. Disadvantages are its higher algorithmic 
complexity, lower numerical stability, the fact that errors are not Galilean invariant, a tendency 
to overmix fluid, and that runtime memory requirements grow with refinement. The most com- 
monly used AMR codes in the galaxy formation community are Hydro- Art 12421 13301 13311 . 
Enzo 13321 13331 . and Ramses [250| (with the Flash code 113341 about to join the fray 113351 ). 

Both SPH and AMR techniques have weaknesses, which are related to the numerical ap- 
proach taken to solve the fluid equations. One advantage of SPH is its pseudo-Lagrangian nature 
which fits very well the needs of cosmological structure formation simulations, where adaptivity 
and a large spatial and dynamical range is required. On the other hand AMR, as a finite volume 
scheme, provides highly accurate results for fluid problems providing, for example, very good 
resolution of shocks, discontinuities and mixing, which are typically harder to resolve very well 



28 



with SPH schemes. A natural way to combine the advantages of SPH and AMR techniques is 
to allow for Moving Meshes in the volume discretization. This idea goes back to the 1990's, 
where moving meshes were first explored in the context of astrophysical applications 13381 13391 . 
Although the idea of having the computational mesh move with the hydrodynamical flow seems 
very natural, its practical implementation turned out to be rather difficult. Approaches relying 
on deformed Cartesian grids lead to problems in handling grid deformation properly in fully 
astrophysical applications 1 3381 13391 . Only recently new moving mesh schemes have been de- 
veloped, which are able to circumvent this problem 113 131 13401 . These new schemes do not use 
coordinate transformations like previous moving mesh codes in cosmology, but instead employ 
an unstructured Voronoi tessellation of the computational domain. The mesh-generating points 
of this tessellation are allowed to move freely, offering significant flexibility for representing the 
geometry of the flow. If the mesh motion is tied to the gas flow, the results are Galilean-invariant 
(like in SPH), while at the same time a high accuracy for shocks and contact discontinuities is 
achieved (like in Eulerian schemes). Furthermore the mesh is free of the distortion problems 
inherent to previous schemes. Using the Arepo code [ 3 1 3 1 , this new computational approach 
has recently been applied successfully to initial large-scale cosmological simulations of galaxy 
formation EHES) (see Fig.g). 

An accurate treatment of the hydrodynamics is necessary, but far from sufficient. In order to 
model the galaxy formation process, simulations must go beyond adiabatic hydrodynamics and 
include gas cooling through radiative energy losses, as well as heating from an externally cal- 
culated meta-galactic UV background 1134 II 13421 . The cooling is typically implemented through 
tabulated cooling functions, which provide externally calculated (using the Cloudy code H3434 - 
345 1, for example) equilibrium cooling rates as a function of density, ionization fraction, temper- 
ature, metallicity, and intensity of the UV background. Some codes follow the non-equilibrium 
abundance of hydrogen and helium species (including molecular hydrogen in some cases) cou- 
pled to gas cooling and heating by solving a chemical network sub-cycled on each hydrodynamic 
time step. The details of how gas cooling is implemented can make a difference in the simula- 
tion's outcome. 

In regions that have cooled and condensed to sufficiently high density, stars will form. This 
process is captured with a sub-grid model, in which a fraction of the available gas is converted 
to "star particles" representing an entire stellar population, which from then on are treated as 
collisionless and evolved in the same fashion as DM particles. The nature of this star formation 
(SF) prescription is another important differentiating aspect of galaxy formation simulations. 
Most commonly, a Schmidt law [ 346 1 is employed, whereby the star formation rate (SFR) is 
proportional to the gas density divided by the local free fall time, resulting in a SFR propor- 
tional to p 1 ! 2 . However, some authors instead prefer a linear SF law corresponding to a fixed 
SF time scale 13301 . especially with simulations that distinguish between atomic and molecu- 
lar gas phases and tie the star formation to the latter 13361 13471 13481 . The SF efficiency is a 
free parameter in principle, but is usually set to a few percent, motivated by observational con- 
straints 13491 . Furthermore, a SF density threshold is often employed to limit SF to especially 
high density regions. Some methods additionally restrict SF to regions with a converging flow 
and a cooling time shorter than the dynamical time [ 350 1 . In order to prevent the formation of 
excessive numbers of star particles, a minimum star particle mass is sometimes enforced. For 
any given model, the parameters are typically tuned to reproduce macroscopic SF scaling laws 
like the Kennicutt-Schmidt relation [351 1. Nevertheless, with so many different SF prescriptions 
and tunable parameters, it is not surprising that there is no unique solution, and that results vary 
greatly, especially in conditions far from where the prescriptions were calibrated. 
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Last, but by no means least, galaxy formation simulations must account for so-called feed- 
back processes, which attempt to capture the injection of mass, momentum, energy, radia- 
tion, and metals (nucleosynthetic products) from massive young stars, evolved asymptotic gi- 
ant branch stars, exploding supernovae (both type la and II), and accreting black holes into the 
surrounding interstellar medium. Feedback appears to be a crucial ingredient in explaining the 
macroscopic properties of individual galaxies [337 352 -3551 as well as galaxy population statis- 
tics |356 -358 ], and may be the key for explaining the small scale challenges CDM faces (see next 
section). At the limited resolution of current simulations, feedback is implemented using heuris- 
tic and often ad-hoc prescriptions (for a review, see 1 359 1 ). Numerical difficulties that must be 
overcome include preventing the injected energy from immediately being lost owing to the high 
densities and cooling rates in star forming regions, accounting for radiation pressure from young 
massive stars, forming and maintaining large scale galactic outflows, and properly accounting 
for the mixing of metals. Even more so than for the SF prescription, results depend sensitively 
on the details of how feedback is implemented [323], and this is probably the greatest source of 
uncertainty in present day galaxy formation simulations. 

For completeness, we also briefly mention here that baryonic physics effects are also com- 
monly accounted for through the use of Semi-Analytic Models (SAM). In this approach gas 
cooling, star formation, feedback from SNe and AGNs are all implemented in a simplified but 
self-consistent manner on top of halo matter merger trees derived from analytical calculations 
or from DM-only simulations. SAMs have been used to explore the connection between dark 
matter structure and galaxies (for a review, see 13601 ) have helped to quantify how uncertainties 
in galaxy formation can propagate to DE studies B167l fl83l. 

3.4.2. Baryonic effects on DM 

Given the uncertainties arising from the choice of hydrodynamic method, the sub-grid physics 
prescriptions, and the large number of adjustable parameters, it is not surprising that there is not 
yet consensus on how baryonic physics alters the distribution of DM in halos. In the following 
we report on some of the results from recent hydrodynamic galaxy formation simulations. We 
caution the reader, however, that in most cases these results should not be viewed as definitive 
answers, and that the conclusions are subject to change with higher resolution and improvements 
in the treatment of sub-grid physics prescriptions. 

Concerning the large scale distribution of DM, it has been shown 13311 13611 1362 1 that the 
inclusion of baryonic physics substantially alters the matter power spectrum on scales (k > 
lh Mpc~' or I > 800) that are relevant for contemporary and future weak lensing galaxy surveys 
aiming to constrain the nature of DE. The weak lensing shear signal has also been shown to 
be affected by baryonic physics, especially when strong AGN feedback in galaxy clusters is 
accounted for 136311 . and this can lead to significant biases (tens of percent) in the inferred DE 
equation of state parameter wq. DM halo mass functions can also be affected, with baryons 
causing a 10% enhancement in the cumulative mass function for > 10 12 /i _1 M halos in one 
study 113 3 1 II . and ~ 30 percent deviations in the number density of 10 H h _1 M halos in another 
136411 . with the sign of the deviation (increase or decrease in halo mass) depending on how 
the baryonic physics was implemented. Note that 10 - 30% differences are larger than the 5% 
statistical uncertainty in DM-only halo mass function 11221 . Gas mass fractions in galaxy clusters 
and halo mass scaling relations of the thermal Sunyaev-Zel'dovich effect depend on the treatment 
of baryons [365 1. Together these results imply that deep galaxy cluster surveys designed to tightly 
constrain cosmological parameters and the nature of DE must account for baryonic effects. 



30 



The abundance and spatial distribution of subhalos inside galactic and cluster scale halos is 
another area likely affected by baryonic physics, and here again even the sign of the effect is 
unknown. Baryonic condensations within subhalos could increase the central density and make 
them more resiliant to tidal disruption. This could lead to an increase in the subhalo abundance 
close to the halo center, as seen in some simulations [366 . 367 1 . On the other hand, if the host halo 
itself has a sizeable stellar disk, then disk shocking H368L as well as interactions with individual 
stars II369L could lead to enhanced subhalo destruction and lower their abundance in the inner 
part of the galaxy. 

From the point of view of DM detection experiments, the most significant concern is the 
possibility that baryonic effects could modify the shape of the DM density profile. This is the 
arena of a long standing debate between advocates of adiabatic contraction increasing the central 
DM density on one side, and proponents of processes removing DM from the halo center on the 
other. If the cooling and condensation of gas proceeds slowly, with baryons gradually sinking to 
the center, then DM will be dragged in adiabatically, leading to a steepening of the DM density 
profile from the NFW slope of ~ -1 to something closer to isothermal -2. This adiabatic 
contraction effect was worked out analytically in 037011 and its description has since been refined 
13711 13721 . It is routinely observed in cosmological galaxy formation simulations 13371 13721 - 
13771 . 

If, on the other hand, baryonic material is rapidly delivered to the center, for example through 
cold flows |378 1, then adiabatic contraction may not operate and other dissipationless processes 
that transfer energy from the baryons to the DM could lower the central DM density [ 379 1 . Ex- 
amples of such cusp-to-core conversion processes are resonant interactions of the DM halo with 
a stellar ba^] 1 382) 43 841 , the decay of a supermassive black hole binary M385I , dynamical friction 
of dense stellar clumps against the smooth background DM halo B3661 13731 13861 [3871 . and re- 
peated and violent oscillations in the central potential due to energy injection from active galactic 
nuclei [388] or supernova-driven galactic outflows [389 -392]. This last process in particular has 
recently attracted much attention, since galaxy formation simulations with very efficient "blast- 
wave" supernova feedback have for the first time resulted in spiral galaxies with realistically low 
bulge-to-disk ratios H3371 [3541. and also exhibit pronounced DM density cores B393I . The re- 
moval of the central DM density cusp through baryonic processes is especially interesting in the 
context of the Missing Satellites and Too Big To Fail problems of CDM (see jgg) 0941 13951 



Lastly, we mention the possibility that baryonic physics effects could lead to a displacement 
between the point of maximum DM density in a halo and its dynamical center [396 , 3971 . as 
recently identified in the Eris simulation 0371 . one of the highest resolution and most realistic 
cosmological hydrodynamic simulations of a Milky Way-like galaxy (see Fig. |4j. In this simu- 
lation, but not in comparable DM-only simulations, the DM offset was ~ 340 pc averaged over 
the last 8 Gyr, typically in the plane of the stellar disk, and aligned to about 30 degrees with the 
orientation of the stellar bar. Such an offset would considerably alter expectations for the indirect 
DM detection signal from the Galactic Center, where Sgr A*, the compact radio source associ- 
ated with a supermassive black hole, is a good marker of the dynamical center of the Galaxy 
II100I1 . Intriguingly, the recently reported 130 GeV gamma-ray line from the Galactic Center 
B105I 4T08I appears to be offset by about 1.5 degrees (about 200 pc projected) from Sgr A*. 



6 Some authors, however, instead find that the formation of a stellar bar actually increases the central DM density 

1380113811 . 
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4. The Next Decade 



Having rextensively reviewed the current state of the art of cosmological DM and galaxy 
formation simulations, we now present our vision for this field for the next ten years. What 
progress should be possible and where should priorities be focused in order to maximize our 
understanding of the DM and DE problems? 

4.1. Dissipationless Dark-Matter-only Simulations 

On the cosmic scale, there is a need for very large numbers of high resolution simulations 
in large box sizes, scanning over cosmological parameters. Future DE surveys (e.g. DES, Big- 
BOSS, LSST) will cover enormous volumes (10000's of square degrees out to z ~ 2) and are 
sensitive to quite low mass galaxies. In order to perform grid-based or Markov Chain Monte 
Carlo estimations of cosmological parameters, their errors, and especially the co-variances of 
their errors, it will be necessary to have highly accurate theoretical predictions of the non-linear 
clustering of matter for a finely sampled scan of cosmological parameters. At the desired per- 
cent level accuracy, such predictions can and must come from cosmological simulations [398 1 
- hundreds to thousands of them [399, 400]. Each simulation must cover volumes comparable 
to the surveys, and have a mass and force resolution sufficient to resolve nonlinear clustering on 
galactic scales. Furthermore, owing to non-linear mode coupling, small scales are affected by 
the sampling variance of the largest modes, necessitating multiple realizations for any given cos- 
mology. In the mildy non-linear regime (k < 0.3/;Mpc _1 ), it may be possible to greatly speed 
up cosmological parameter scans by utilizing rescaling algorithms |268| or second order La- 
grangian Perturbation Theory in conjunction with time- and scale-dependent transfer functions 
that provide good approximations to particle trajectories and can be derived from cheap N-body 
simulations 1140111 . 

Beyond the standard ACDM simulations, there is considerable interest in exploring simula- 
tions with alternative gravity laws, such as f(R) theories, that explain cosmic acceleration without 
DE. In addition to the N-body gravity treatment, such simulations must also solve for the non- 
linear dynamics of the Chameleon scalar field that screens the gravity law modifications from 
small scales. In many cases, this is accomplished with adaptive multi-grid relaxation techniques, 
quite similar to how gravity is treated in AMR hydrodynamics codes. Promising initial progress 
is already being made in this area 1402 14051. 

On cluster and galactic scales, efforts will be focused on the substructure problem, and will 
attack the problem from three directions: 

i) Higher resolution 

Given that the CDM substructure hierarchy extends for 10 - 15 orders of magnitude below 
current resolution limits, direct numerical simulations will realistically not be able to resolve 
the full hierarchy even in the intermediate future. Nevertheless, pushing on resolution in 
zoom-in simulations of individual halos (as in the Via Lactea and Aquarius projects) is a 
worthwhile goal, for the following reasons: 

• Stars in the newly discovered ultra-faint dwarf galaxies of the Milky Way are typically 
only found out to < 100 pc from the dwarf's center. In order to compare stellar 
kinematics in these systems with predictions of central densities and profile slopes 
from CDM simulations, it will be necessary to push to a mass and force resolution 
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of ~ 100M o and ~ 10 pc. Note that baryonic physics modifications are expected to 
be less important in the low mass (< 10 9 M Q ) host halos of ultra-faints 1139011 , and so 
DM-only simulations should still provide valuable predictions. 

The DM annihilation boost factors from substructure depend sensitively on the prop- 



erties of the subhalo population below current simulation's resolution limit (t 2.2.2 ix). 
Since the power spectrum of DM density fluctuations is not truly scale invariant, one 
might expect quantitative changes in subhalo properties at lower masses. Being able to 
resolve the abundance, spatial distribution, and density profiles of < 10 5 M Q subhalos 
would help to clarify the relevance of substructure boost factors. 
• The phase-space structure in the local neighborhood is only beginning to be resolved 
by current simulations. Higher resolution simulations will provide a somewhat more 
fine-grained view, and will additionally resolve more tidal streams from disrupted sub- 
halos. This will allow a better assessment of the importance of velocity substructure 
for the intepretation of DM direct detection results. A caveat is that baryonic affects 
are likely to be very important. 

The "Silver River" simulation is an example of an on-going effort in this class. With a 
particle mass of 100 M and a force softening of 27 pc, this simulation is pushing the frontier 
in Galactic zoom-in simulations. The run has progressed to z ~ 5 and will be completed 
to z = in the next 1-2 years, pending computational resource support from national 
supercomputing centers (30 - 40 million core-hours required). 

ii) Different host halos 

The total number of ultra-high resolution (> 10 9 particles) simulations of individual ha- 
los is still quite small: one cluster simulation (Phoenix) and three galactic scale ones (Via 
Lactea II, Aquarius, and GHalo), and all four of these simulations have been run with 
only two codes, Pkdgrav2 and Gadget. The six lower resolution Aquarius simulations and 
Millennium-II notwithstanding, the question of cosmic variance, of how much halo-to-halo 
scatter there is in the substructure population, is not fully settled: see for example [406 1, 
who find a considerably larger halo-to-halo scatter in the subhalo abundance than what is 
seen in Aquarius and Millenium-II. Trends with host halo mass and accretion history are 
equally important and in need of further study. These questions are of particular relevance 
to the Too Big To Fail problem (jj |3.2| i, where even a factor of a few variation would go a 
long way towards a resolution B280I . Lastly, current simulations don't properly account for 
the immediate environment of the Milky Way, because they typically don't have a nearby 
massive M31-like companion. 

The Rhapsody Cluster Resimulation Project [407 1 is an early effort along these lines, con- 
sisting of 96 zoom-in simulations of cluster scale halos with ~ 10 7 particles per halo. Sim- 
ulation suites of hundreds of Milky Way-like halos with varying masses and accretion his- 
tories are also currently being pursued. Capturing more of the Local Group structure can 
be achieved either through constrained realization simulations [408] or by picking suitable 
halo configurations from an ensemble of lower resolution simulations (S. Garrison-Kimmel 
et al., work in progress). 

iii) Alternate DM models 

Going beyond the cold and collisionless DM paradigm, sharpening the predictions that 
WDM or SIDM models make will be another promising avenue of exploration. In models 
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with a cutoff in the density power spectrum it will be interesting to study how halo forma- 
tion is suppressed right around the free streaming scale. This is numerically challenging, 
due to the extremely high resolution required to push the artificial fragmentation to small 
enough scales. If this challenge can be overcome, it may be feasible to conduct ultra-high 
resolution Milky Way-scale simulations (like Via Lactea II or Aquarius A-l) with one of the 
alternative DM models. This would be very useful for comparing with stellar kinematic data 
in MW dwarf satellite galaxies and with data from future strong lensing surveys of flux ratio 
anomalies. It may also be of interest to simulate not generic WDM or SIDM, but specific 
well-motivated particle physics models, for example the Atomic DM model |13| with its 
baryonic-like small scale power spectrum oscillations. 



4.2. Simulations Including Baryons Physics 

We anticipate that improvements in the treatment of baryonic physics processes will be a ma- 
jor focus in cosmological simulations over the next decade and beyond. As detailed in 5 3.4 there 
are numerous conceptual and technical challenges in properly implementing baryonic physics. 
Yet these are important problems to tackle, since they have potentially far reaching consequences 
for our understanding of the distribution of DM in and around galaxies, groups, and clusters, and 
directly affect expectations for indirect and direct detection experiments. 

On cosmic scales, it will be crucial for numerical simulations to quantify how baryonic 
physics modifies the matter power spectrum 113631 , how it affects the bias between galaxies 
and DM halos, and how these effects impact cosmological parameter estimations from upcom- 
ing surveys. For individual cluster simulations, it is imperative to establish how reliable X-ray 
and Sunyaev-Zeldovich measurements of halo masses are |409|, and how well cluster mass- 
observable relations can be calibrated. In addition to radiative cooling, star formation, and stellar 
feedback, such simulations must account for magnetic fields and anisotropic conduction, non- 
thermal pressure support from turbulence and cosmic rays, as well as AGN feedback. 

On galactic scales, one of the most urgent questions that numerical simulation must strive 
to clarify is under what conditions adiabatic contraction steepens the central DM density pro- 
file, and whether this effect can be overcome by feedback processes that may redistribute large 
amounts of DM and result in shallower or even cored density profiles. If both processes occur in 
nature, which one dominates, and how does the answer depend on halo mass, environment, and 
cosmic time? It seems clear from past work that simulations with cooling, but little or inefficient 
supernova (SN) feedback, exhibit a substantial amount of adiabatic contraction of the DM halo. 
On the other hand, these simulations typically also suffer from a baryonic over-cooling problem, 
and produce too many stars, that are too centrally concentrated. Some form of stellar feedback 
thus appears to be a necessity. Whether the resulting regulation of star formation is accompanied 
by a removal of substantial amounts of DM from the central regions is an open question. 

When efficient SN feedback is imposed "by hand", either by artificially turning off gas cool- 
ing, or by employing unphysically high star formation or SN energy injection efficiencies, or by 
hydrodynamicaily decoupling wind particles from the surrounding gas, substantial redistribution 
of DM has been observed in a variety of simulations. But how realistically are these ad-hoc 
implementations capturing the actual physical processes occuring in star forming regions inside 
giant molecular clouds? Is the resulting DM cusp flattening a robust outcome? The answers await 
more sophisticated treatments of star formation and feedback processes. Promising directions of 
future investigations along these lines include, but are not limited to: 
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• Suppressing star formation in low mass and low metallicity systems by using Fb-regulated 
star formation prescriptions rather than (or in conjunction with) SN feedback 1 3361 13471 : 

• Accounting for radiation pressure from young massive stars [410|, which imparts momen- 
tum into the surrounding gas and can increase the efficiency of subsequent SN feedback 
without resorting to artificial enhancement; 

• Modeling non-thermal support provided by unresolved turbulence, magnetic fields, and 
cosmic ray propagation 1141 1114121 . 

In general more sophisticated treatments of star formation and feedback physics will require 
resolving the hydrodynamic component of the simulations with parsec scale resolution (the DM 
can be treated with coarser resolution for these purposes). This is at least a factor of ten higher 
than what is currently achievable in cosmological zoom-in simulations. Full-box simulations 
with this resolution will be limited to the high redshift domain, and thus it will be challenging to 
obtain large samples of realistically simulated galaxies for comparison with local observations. 

It is commonly hoped that it will eventually be possible to include realistic star formation and 
feedback prescriptions in low resolution, large volume, full-box simulations by calibrating them 
to the results from much higher resolution zoom-in, or even isolated, simulations, in which a 
more accurate treatment of the relevant physics is possible. However, it is yet to be demonstrated 
that this approach will be feasible. Recent algorithmic developments in creating initial condi- 
tions (by generating white noise fields hierarchically in real space [220, 266 1) are now enabling 
simulations of unprecedented dynamic range. The MUSIC code [ 266 1 has demonstrated high ac- 
curacy in reproducing in nested simulations the same detailed structures as in corresponding high 
uniform resolution full-box runs. This provides the basis for the Cosmic Renaissance (CORE) 
Project (T. Abel et al., work in progress), which aims to conduct a large set of self-consistent 
nested simulations over a wide range of scales, from extreme zoom-ins on the formation of the 
first stars in a volume the size of the observable Universe to increasingly less focused and more 
coarsely resolved simulations of the formation of galaxies, clusters, and the large scale structure 
of the Universe. 

Lastly, it is also necessary to understand how these numerical results depend on which hy- 
drodynamic technique and indeed which code was used. The Aquila Code Comparison Project 
B3231 was an important step in this direction, but owing to its approach of allowing simultaneous 
variations in code and feedback physics implementations, it was difficult to disentangle which 
factor the observed differences can be attributed to. Nevertheless, it clearly demonstrated that 
the current generation of cosmological hydrodynamic simulations are not yet able to uniquely 
predict the properties of the baryonic component of a galaxy forming in a fully specified dark 
matter accretion history. The recently initiated Santa Cruz High-resolution Galaxy Simulation 
Comparison ProjecQis a similar effort, focusing on < 100 pc resolution galaxy formation simu- 
lations with a wide range of current state of the art codes representing SPH, AMR, and Moving 
Mesh techniques. 

4.3. Computational Trends and Algorithmic Advances 
4.3.1. Processing 

Over the last decade there have been significant advances in High Performance Computing 
(HPC) available to astrophysical research. The top 20 machines in the June 2012 TOP500 lis{^] 



7 https://sites. google.com/site/santacruzcomparisonproject/ 
s http://www.top500.org/list/2012/06/100 
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have all exceeded the petaflops barrier, and utilize upwards of 100,000 cores to do so (Sequoia, 
the current No.l has 1,572,864 cores). Supercomputers are expected to reach the exaflop scale 
during the next decade. This will enable us to carry out much more detailed modeling of the 
problems described in this manuscript, but it also poses serious challenges for the development 
of numerical codes and algorithms capable of fully exploiting new technologies. 

i) Scalability and performance 

One characteristic of current and future supercomputers is the very large number of comput- 
ing cores. The most obvious problem concerns how to distribute the computational load over 
an increasingly large computational domain, while at the same time keeping communication 
time to a minimum. Preserving adequate balance in the CPU and load requirements will re- 
quire development and improvements in all algorithms present in N-body calculations, in 
order to achieve scalability to millions of compute tasks. 

Another interesting aspect of supercomputing in the next decade will be the architecture of 
their Central Processing Units (CPU). Modern compute nodes contain multiple processing 
cores, each of which have their own sets of instructions and cache, but have access to com- 
mon memory (RAM). In the near future we can expect nodes with hundreds of such cores, a 
feature that needs to be exploited by simulation codes via mixed parallelization schemes, in 
which the commonly used Message Passing Interface (MPI) model for distributed memory 
is combined with shared memory parallelism, via e.g. OpenMP or Posix threads. Some ex- 
isting N-body solvers already take advantage of SIMD (Single Instruction Multiple Data), 
which allows parallel vector calculations in the CPU. Speed-ups can be factors of a few, 
however exploiting SIMD requires machine-dependent code. 

In recent years, Graphical Processing Units (GPU) have become a competitive alternative 
to CPUs for intensive numerical tasks. Initially developed for rapid and highly parallelized 
manipulation of computer graphics, they have also found wide use for general-purpose com- 
putation. For certain algorithms able to take advantage of data-parallelism, GPUs can result 
in very significant speed-ups of factors of ten or more. HPC is already beginning to embrace 
this trend, with 3 of the top 10 supercomputers (Tianhe-IA, Jaguar, and Nebulae) taking ad- 
vantage of GPUs. A downside of using GPUs is the comparatively small amount of memory 
available on board, requiring a lot of slow data transfer from CPU to GPU and back. Nev- 
ertheless, there have been efforts to implement widely used algorithms to take advantage of 
GPUs, for instance Fast Fourier Transforms, kd-Trees, as well as visualization and analysis 
tools. Custom code modifications can be performed using extensions to the C programming 
language, like CUDA (specific to NVIDIA GPUs) and OpenCL (multiple architectures). 
As graphics cards improve in reliability and performance, cosmological N-body codes will 
likely take advantage of this technology (see e.g. [413-4-15 1), most likely in programming 
models that employ hybrid algorithms for both CPUs and GPUs. 

Perhaps a point of convergence between CPUs and GPUs will be in hardware such as the 
Intel MIC cards (Many Integrated Core). This product, resulting from Intel's Larrabee re- 
search project, is essentially a GPU-CPU hybrid and functions as a co-processor for HPC. 
It is expected to become available by the end of 2012 and could become an alternative to 
GPUs with a simpler programming model. 

ii) Big Data 
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As we have seen, state of the art simulations already generate petabytes of data. Handling 
this data is becoming an increasingly difficult problem and we can only anticipate that it will 
get worse during the next decade. One limiting aspect is the disk I/O speed; even with paral- 
lel distributed filesystems (e.g. Lustre) it does not scale with the size of the supercomputer, 
and thus disk I/O can take a considerable fraction of the runtime of a simulation. Another 
aspect is how to manage and store extremely large datasets, especially with collaborations 
commonly spread across the world and the considerable cost of dedicated storage hardware. 
Simply migrating the data from large supercomputing centers to smaller local facilities can 
be difficult and time consuming. 

These are serious challenges for computational cosmologists, and strategies will need to 
be developed to cope with these issues. One natural development already in practice is to 
merge analysis and visualization software with the simulation code so as to reduce as much 
as possible the output of data during runtime and the handling afterwards^] For large DM 
simulations this may mean identifying and extracting (sub)halo information, constructing 
merger trees and lightcones on-the-fly, and storing only the resulting reduced data. However, 
it is never possible to anticipate all the uses of a simulation, and this motivates saving at least 
some full outputs. Novel data compression techniques are a promising way forward. The 
idea is that data analysis often does not require the same precision as run-time computation. 
Substantial compression of the output can the be achieved either by neglecting irrelevant 
bytes or by spatial averaging (a compression approach widely used in image and video 
displays). For instance, the data can be stored in a hierarchical tree fashion, where the 
initial bytes provide a coarse representation of the dataset, and subsequent bytes provide 
refinement in specific quantities and/or spatial regions. These would alleviate access and 
reading time for postprocessing, at a cost of giving up the ability to restart simulations and 
with file formats very specific to a given simulation. 

The distribution of data worlwide is a serious but important problem, as serving data prod- 
ucts is necessary and desirable for scientific progress. In the last few years, an increasingly 
popular option has been the developments of internet databases providing reduced data prod- 
ucts, for example from the Millennium and MultiDark simulation database^] Another al- 
ternative is to store the data itself in a distributed fashion on the internet in the "cloud" and 
using Hadoop for distributed data analysis. 

4.3.2. Novel Approaches 
i) New Gravity Solvers 

Despite well-known problems with the N-body treatment of cosmological structure forma- 
tion (e.g. two body relaxation and artificial fragmentation), there has been little progress in 
how to solve the collisionless Boltzmann equation with self-gravity (Poisson-Vlasov equa- 
tions). As the accuracy demanded of cosmological simulations increases, it is likely that 
the short-comings of current N-body schemes will eventually become limiting. This moti- 
vates improvements in existing treatments, such as the introduction of adaptive gravitational 
softening 14171 14181 or improved density estimators based on following the distortion of a 
tetrahedral tessellation of the initial density field 114 191 14201 . as well as radical departures 



'This is also one of the stated development goals of the yt Project 14161 . http://yt-project.org/ 
°http://www.mpa-garching. mpg.de/millennium/ and http://www.multidark.org/MultiDark/ 
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from the Monte Carlo N-body approach, such as full six-dimensional Vlasov solvers [421 1 
or employing kinetic theory [422 1. 

ii) New Programming Languages - the need for parallelism 

Parallelism and HPC become increasingly important also for industry and business applica- 
tions, where server farms and computing cluster with tens of thousands of cores are now used 
for large-scale data mining and reduction problems. Furthermore, since the break-down of 
Moore's law for the raw core processor performance, developers in all fields, ranging from 
mobile devices to high performance cluster, need to start thinking parallel. It is therefore 
not surprising that new programming languages, paradigms, and models are now emerging, 
which have parallelism more naturally built in. 

One example is Unified Parallel C (UPC) as an extension of the ISO C 99 programming 
language designed for high-performance computing on large-scale parallel systems. This 
includes architectures with common global address spaces (SMP and NUMA) and also dis- 
tributed memory systems, which are more common for cosmological applications. The main 
idea of UPC is to present the developer with a single shared, partioned address space, and a 
Single Instruction Multiple Data (SIMD) programming model. UPC is a specific implemen- 
tation of the more general class of PGAS (partitioned global address space) programming 
languages. Other languages which follow the PGAS principles are Co-array Fortran, Ti- 
tanium, Fortress, Chapel, X10 and Global Arrays. PGAS is based on two main concepts: 
multiple execution contexts with separate address spaces and access of memory locations on 
one execution instance by other execution instances. It is very likely, that in the near future 
simulation codes will start using some of these language. 

Besides the actual simulation software, data processing software also becomes increas- 
ingly complex and requires efficient parallelization. Large parts of available post-processing 
pipelines are programmed in general-purpose, interpreted high-level programming languages 
like Python, which is slowly replacing older IDL implementations. Both languages were es- 
sentially designed without a particular focus on parallelism, though it is possible to use 
them for parallel processing. But new languages are also arising in this field, which will 
likely replace Python or IDL as the de-facto standards in the near future. For example, 
Julia is new a high-level, high-performance dynamic programming language for technical 
computing, which provides a sophisticated compiler, distributed parallel execution, and an 
extensive mathematical function library. Altough it is a LLVM-based just-in-time (JIT) 
compiler based language it reaches nearly C/C++ performance due to heavy optimization 
outperforming IDL, MATLAB, R, Python, Octave by a large factor. Most importently it 
has parallelism naturally build in, which makes it suitable for large-scale data mining and 
reduction on huge compute clusters. 

UPC, as a PGAS parallel programming language, and Julia, as a high-performance dynamic 
programming language, are just two examples and it is very likely that the programming 
language landscape will move significantly more towards HPC and parallel computing in 
the very near future mainly driven by industry needs. 

5. Conclusions 

Two of the greatest mysteries of contemporary physics are the nature of dark matter and of 
dark energy. As we have seen over and again in this work, experimental and observational efforts 
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to get at the answers to these questions are intimately tied to predictions from cosmological DM 
simulations. The simulation field, however, is rapidly evolving, as computational resources con- 
tinue to grow, enabling larger, more complicated, and more realistic simulations to be performed. 
As a result, it has become more difficult for physicists not directly involved with simulations to 
follow the progress of the field, and to stay up to date with the evolving implications for stud- 
ies of DM and DE. Motivated by this realization, we have in this review attempted to provide 
an overview of the current state of the field of cosmological DM simulations, with a particular 
emphasis on the connection to DM detection experiments and observational probes of DE. We 
have summarized the successes and accomplishments of the current generation of simulations, 
but also detailed their problems and short-comings, as well as the challenges faced by the next 
generation of simulations in the coming decade. 

For the DE problem, the greatest need appears to be a tremendous increase in the number 
of very high resolution, very large volume, DM-only simulations, scanning over cosmological 
parameters in order to allow brute-force Monte-Carlo estimates of the parameter errors and their 
co-variances for future DE surveys. At the same time, cosmic scale simulations will continue to 
investigate the effects of baryonic physics, non-Gaussian initial conditions, and modified laws of 
gravity. 

For the DM problem, on the other hand, it seems that the DM-only approach by itself is 
nearing the end of its usefulness - baryonic physics effects are simply too important to neglect on 
galactic scales. There is some room for improvement in the DM-only approach, in particular for 
the internal structure of subhalos, for exploring cosmic variance and host halo mass dependency, 
and for departures from the cold and collisionless DM assumption. However, the main driver of 
the field must be studies of baryonic effects on the distribution of DM in halos, since these have 
the potential to profoundly alter expectations for DM detection. 

The next decade promises to be filled with exciting challenges and the potential for great 
discoveries. It is safe to say that the field of Numerical Simulations of the Dark Universe will 
not run out of things to do. 
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